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ABSTRACT 





A shell formulation was developed from a three-dimensional solid. The shell 
element has four corner nodes at which there are three displacements and three rotations as 
nodal degrees of freedom, and includes both transverse shear and transverse normal 
deformations. The element utilizes reduced integration along the in-plane axes and full 
integration along the transverse axis. The formulation incorporates the Gurson constitutive 
model for void growth and plastic deformation. An algorithm for stable solutions of the 


nonlinear constitutive equations is also developed. Hourglass mode control is provided by 





| adding a small fraction of internal force determined through full integration along the in- 
! plane axes and reduced integration along the transverse axis. Implementation into a research 
finite element program is discussed. Numerical examples are provided to verify the accuracy 
of the new element and to show the importance of the transverse normal stress, void effects 


on plastic strain, and the necessity of applying a drilling moment 
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I. INTRODUCTION 


Since shell structures are efficient load-carrying structural members, they have been 
dominant in most structural applications. However, the finite element formulation of shells 
poses some difficulty. As a result, extensive study has been devoted to developing better 
shell elements. Because of the abundance of papers on this subject, no attempt is made here 
to summarize all of them. Some of the related past work is given in references [1-12]. 

Void growth and nucleation can have a significant effect on plastic flow [13]. Since 
voids act as stress concentrators, the overall effect is to reduce the sii under plastic flow, 
and increase the plastic strain [14-16]. The model proposed by Gurson [17] appears to have 
been adopted as the standard for incorporating void growth and nucleation effects into a 
numerical solution of elasto-plastic problems. Previous work has been done on improving 
the efficiency a accuracy of plasticity computations, applying plasticity to plate/shell 
elements, and incorporating void growth and nucleation effects in solid elements [1 8-20]. 

The element presented in this paper incorporates a modified version of the algorithm for 
three-dimensional solids proposed by Aravas [21] into a shell element. This shell element 
assumes a modified plane-stress condition, and utilizes both the hydrostatic pressure and the 
deviatric stress. Due to the importance of the hydrostatic pressure on void growth and 
sail this element includes the transverse normal ress 

The Radios proposed by Aravas is not unconditionally stable when applied to this 
| modified plane stress condition [21]. This paper also presents modifications to the original 
algorithm, which greatly enhance the stability of the solution, at the cost of a slight decrease 
in computational efficiency. This algorithm also allows for more eompleaied work- 


1 





hardening profiles than the typical exponent law (o = Ke"). Full modeling of the entire 
stress-strain curve is accomplished by a piece-wise linear approximation. While this method 
is not particularly useful for analytic approaches, it appears to be helpful in strictly numerical 
solutions. The stress-strain relationship may be taken directly off the results of a standard 
tensile test. 

For thick shell applications, the transverse normal stress and strain cannot be 
neglected. The transverse normal stress and — affect both the hydrostatic pressure and 
‘the deviatoric stress, which in turn affect plastic deformation and void growth and 
nucleation. The element saciaadal here extends the typical use of the drilling degree of 
freedom (DOF), as outlined in Hughes and Brezzi [22], to incorporate the effects of 
transverse normal strain. The drilling DOF is important in transmitting stress and stein to 
elements across sharp bends and curves, and is therefore already included in a variety of shell 
elements [2,4,6,7,10]. In addition to this traditional usage, the present element utilizes the 
drilling DOF to compute the transverse normal deformation. The importance of including 
the transverse normal deformation is outlined in Essenburg [23]. 

The present study formulates a shell element for transient analysis. The element can 
have elasto-plastic deformation with void growth and nucleation. Gurson's void model is 
used as a basis for the void constitutive model. The shell element includes both transverse 
shear deformation and the transverse normal deformation for thick shell applications. The 


drilling degree of freedom is used for computing the deformation through the thickness of 








a thick shell. An algorithm for stable solutions of the nonlinear constitutive equations is also 
developed. 

Some example problems are presented to evaluate the formulation and to investigate 
the effects of the transverse ee strain for thick shells in association with elasto-plastic 
deformation, including void efiects. The sections that follow are: the formulation of the shell 
element, discussion of the elasto-plastic constitutive equation including void ancieation and 
growth, the stable solution algorithm including spurious (hourglass) mode control, numerical 


examples and discussion of results, and conclusions. 




















I. FINITE ELEMENT FORMULATION 


A. GEOMETRY 


A point in a shell structure can be expressed by a vector sum of two ecto The first 
vector is a position vector from the origin of the global coordinate system to a point ona 
reference surface of the shell element. The second vector is a position vector from this 
reference surface to the point under consideration. The surface that spans the center of the 
transverse axis is used as the reference surface in this formulation, although any surface 
would suffice. The first vector terminates at the reference surface directly below the point 
in question. The second vector is then the normal from the reference surface that intersects 


the desired point. Figure 1 shows this relationship. Two shape functions are used to describe 
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Figure 1. Element Cross Section. 
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a position in the element; N~ is the two dimensional shape function in the E-n plane, and 
H “is the one dimensional shape function along the C axis, where (£,n,C) describes a point 
in the natural coordinate system. A generic point in the shell may now be described in terms 


of the position vectors of the nodes and the shape functions: 
uN C= ON (En) +N EM HOG (= 12,3) (1) 
k=] k=} 


where x; is the position vector of node k in the reference surface; )’%, is the unit vector at 


the node k; and n is the number of nodes per element. In the present formulation, a four-node 


shell element is considered. The unit vector ’%, is defined as: 


(x* \ — (x! ) aaa 


3i Ix? \? — (x! aaa 


(2) 








where top and bottom indicate the top and bottom surfaces of the shell, and || || denotes the 


Euclidean norm. The one-dimensional shape function H Kis expressed as: 





(3) 





H"(¢)= F (l+¢)0-¢) -7d ~¢)\(1 +2)| Int)” — (xf ye" 


in which ¢ indicates the location of the reference surface and varied from -1 to 1( ¢ =0 


denotes the mid-surface). The two-dimensional shape function Nk is expressed as: 


N'=4(1-E)(1-n) 


N?=4(1+E)(1-n) . 
NP=F(1+E)(1+0) 
N4=4(1-E)(1+n) 











B. DISPLACEMENT 


The displacement field in a shell can be written as: 
ui(E..o) => N* (En) ul + ON (En) HS) AVG AVS +VEOS) G=12,3) (5) 
k=] k=l 


in which u; is the displacement along the x; axis, 1 is the nodal displacement at the node 
k, and unit vectors Vj; andj/5, lie along the reference surface. V1;,/5, and/’%, are mutually © 
perpendicular. 6} . g; and @3 are rotational degrees of freedom along the unit vectors 


ViisV>; andy, respectively. The right-hand rule is assumed for the positive direction of 


each rotation. Figure 2 illustrates the relationship among these vectors. 





Figure 2. Displacement Vector Orientation. 








g; and 63 are the bending rotations, while g§ is the drilling rotational degree of 


freedom. The role of g@} can be clarified by considering a flat plate parallel to the x,-x, 


plane. Now Eq. (5) can be rewritten for the transverse displacement as: 


u;= > N“uy +> N* HOS (6) 
k=! 


k=1 


Equation (6) demonstrates that the transverse deformation varies through the plate thickness 
(i.e. along the ¢ axis) with 6. In this way, the transverse normal deformation is included 


in this formulation, along with the transverse shear deformations. 


C. COORDINATE TRANSFORMATION 


Combining the three unit direction vectors into matrix [T,] provides the rotation 
transformation matrix, or matrix of direction cosines, as shown in Eq. (7). [T,}° is used to 
transform the nodal degrees of freedom from the global coordinate system to local 
coordinates, as shown in Eq. (8), where k is the node number. The components of {d} at the 


four nodes of an element are shown in Eq. (9). Once the internal force vector is generated 


in local coordinates, [T,] is used to transform the local vectors back into global coordinates. 


This procedure will be discussed later. 


Ir, | Vi V2 Vs las) (7) 








Il 0 0 0 
0 10 0 
00171 0 
{d‘hen*|9 9 0 {dota fen K=1,2,3,4) (8) 
900 pf 
0 0 0 
{dd={d'aaay (9) 


The strain transformation matrix, [T], is used to transform calculated strain from the 
global coordinate system to local coordinates. Transforming the resulting stress from local 
coordinates to the global coordinate system would normally require using [T]", but since [T] 
is orthogonal, [T]’ = [T]", where [T]? is the transpose of [T]. This property negates the 
requirement to invert a six by six matrix. For a detailed derivation of these transformations, . 
refer to Cook [24]. The strain transformation mate is explicitly defined in Eq. (10), where 


V;; is the cosine of the direction vector V; in the x; direction. 


Vi Viz Viz Vir Viz Viz Vis Vii Vis 
Vi Vin V23 Va V2 V2 V3 Vo1 Vas 
2 : : Vn V Vio Vs Va V 
ir]= V31 V32 V3s 31 V32 32 V33 31 V33 (10) 


2ViVa 2Vi2V2 2Vi3sV2 VuVe2tVaVi2 V2VatV2Vis Vis Vat Vo Vir 
2VuV31 2V2V32 2V23V33 Var V32 + 'V31 V22 22 V33 + 'V32 Vos. V23 V31 + V3 V2 
2VuV31 2Vi2V32 2Vi3V33 Vir V32t'V31 Viz Vi2 Vas + V32 Vis. Vis Vai + V3 Vin 





D. STRAIN DISPLACEMENT RELATION 


The six components of the strain tensor are computed from Eq. (5) by taking its 


derivative with respect to the xj sce Tianaiie een anemone elements 
e}= [B] ia} (11) 
fe}={ en ene ¥n%aVis} (12) 

where | 
‘[B)=[[B'][B’ |B B*]| (13) 


The detailed expression for [B‘] is: 











nw 0 - git, og Vi giVis 
0 0 -85Vb 2Vi, gV 52 
0 0 - 2553 g3V43 g3V is 
[B* |= | (14) 
a a O - g5 Vir- 2; Vi» g Vint 2; Vin gVitg, V52 


0 we we _gtyt-gtvt:s ViteVs giVote V's 





a - 83V51- BV g,Viitg,Vis gViit 2, Vis 


in which 
k k 
gt = ONT pty yh OH (15) 
Ox; Ox; 
The vector {d* } is defined as: 
{at }={ ububa§ OF OF OF | (16) 
10 














where uj is the displacement along the x; direction at node k, and © is the rotational 
displacement about the global x; axis at node k. The matrix [B‘] must be calculated for each 


integration point. 


E. JACOBIAN MATRIX 


Computing the derivatives 2 ~ and 24 ~ requires the Jacobian matrix, defined as 


Xre X2E X3¢e 
[J]= Xin Xan Xz | (17) 


X1G X26 X3¢ 


where 
Ox; GON ON" 
X,.=—=) —-x, + H'V3; G@=1,2,3) (18) 
"Oe ty Oe - ae 
Ox; oe 7 Fe) 
in xi ays a ye N" H’V3, G@=1,2,3) (19) 


On ja, ON yl ON 





ng =a wy, 12.9 (20) 
k=] 


[R] is defined as the inverse of the Jacobian matrix, [J]’. Then the required partial 


derivatives are defined as: 

















0 N* ON* ON. | | 
= Ki + R; =1,2,3 21 
ax, = Ri az Ri 5 (i ) (21) 
6 Ht OH* : | 
a i =1,2,3 $9) 
a Ris ac (i ) (22) 
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F. STRESS-STRAIN RELATIONSHIP 
The strain calculated in Eq. (11) is in the global coordinate, and is transformed to a 
local coordinate system using 
{ eiocar $= [ T] egiota } (23) 
where [7] is defined in Eq. (10). Stress is calculated from the strain a the local coordinate 


system using the plane-strain assumption for the in-plane stress components: 





E | =G 
Ox 7 s(e:+V ey) ee 
“VY 
=KG (24) 
oy F s(ve.te,) a Vie 
“V 
Bl ta=KGy,, 


where K is the shear correction factor, E is the elastic modulus, G is the shear modulus, and 

v is Poisson’s ratio. The resulting stresses are in the local coordinate system and are 
converted to the global coordinate system with 

io global }— [TI : { © local } (25) 

{o}=loxoyo2%q Veta} 6 


where [T] is from Eq. (10). 


G. DRILLING MOMENTS 
Let u3 be the transverse deflection as defined in Eq. (6). The work done by a 
pressure loading, p, on an element is: 


W =| u, pda (27) 
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where A® is the element area. Substituting Eq. (6) into Eq. (27), the work is now: 


W ={u,}' [ [IN]p a4+{0,} [ [NH] a (28) 
The first term gives the conventional forces at each node, while the second term yields the 
new nodal load, which will be called the Drilling Moment (DM). For example, if p is a 
concentrated force P at node n, then the drilling moment associated with 6; becomes 
¥%¢ tP , where tis the shell thickness and the mid-plane is the reference plane. 
When P is applied at the top plane, ¢ equals one, and the drilling moment is WtP. 
If P is applied at the bottom plane (still with a positive loading direction), the drilling 
moment is — 4¢P. That is, the load results in compression or tension in the transverse 
normal stress depending on whether the loading is applied to the top or bottom surface of the 
shell. The transverse normal stress is assumed constant through the shell thickness for elastic 
deformation. Ghee as plastic deformation progresses, the transverse normal stress 
becomes non-uniform through the thickness in order to satisfy the yield function. 
The drilling moment direction coincides with the direction of drilling rotation, but 
affects the transverse normal stress. The drilling moment is used as a mathematical 


convenience, and thus, lacks a physical interpretation. 
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H. INTERNAL FORCE, MASS AND THEIR ASSEMBLY 


The stress resulting from Eq. (25) is then converted to an internal force vector and 


summed over all integration points using 


{ ta J if BI’ {o}dédndg=S YY BE fo}wiwewlT] 29) 


isl j=l k=l 


ro Coomera 


where | J | is the determinant of the Jacobian matrix, mx, ny, and nz are the number of 
integration points in the ¢, 7, and ¢ directions, respectively, and w;, w;, and w%, are the Gauss 
weights related to those integration points. The resulting force vector, | Fm {> Must have its 


components transformed back to the global coordinate system using: 


lO 0 0 0 
0 1 0 0 0 
: 00 7] 
{Fv hen=lo 9 9 { tke ben (K=12,3,4) (30) 
000 {rl 
0 0 0 


A lumped mass method is used for the element mass matrix. The matrix for each 


element is diagonal, with equal diagonal elements: 


m OO OO. O 
0 m oO... O 
[M.l=| 0 Om. 0 (31) 


where 
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nx ny nz 
Ww wel J 
[3 2 oy wim me J. (32) 


Hie = 
with 7 = 4 in this four-node element. Using a diagonal mass matrix greatly simplifies time 
integration, since inverting it is a trivial task requiring little computation time. The internal 
force vector and element mass vector (each 24 by 1) are then assembled into the 


corresponding system vectors. 


I. EXPLICIT TIME INTEGRATION 


The use of internal force vectors and explicit time integration negates the need to 
explicitly form the system stiffness matrices. The acceleration vector is computed from 


fo }=[MI' ({ Fe} - {Fe} ) 33) 


where 1U } is the system acceleration vector, [M] is the system mass matrix, ; Fext } is the 
system external force vector, nA superscript ¢ denotes the time step. Of course, the system 
mass matrix in Eq. (33) is simply symbolic, since a system mass vector, formed from the 
diagonal of the mass matrix, is used in actual computation. Velocity and displacement are 


then found using 


(34) 
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lil. DAMAGE CONSTITUTIVE EQUATIONS 


A. GURSON'S VOID MODEL 


Yielding and plastic deformation in the element follows the model proposed by 


Gurson for symmetric deformations around a spherical void [17]. The yielding condition is 


; | 
r-(4| + 24,@cosh{ 322) - (7 + q,7 )=0 (35) 
O00 2 Oo 


where @is the current porosity, p is the hydrostatic stress, q is the effective stress, and o, 
is the current yield stress. This model assumes equivalent yield stress in both tension and 
compression. The constants g,, g>, and q; were introduced by Tvergaard in order to provide 
a better match with numerical studies [16]. Aravas provides a detailed explanation of 
implementing this model in a static finite element algorithm for three-dimensional solid 
elements [21]. The procedure used here is essentially the same, with some modifications due 
to the different element formulation. The stress is then transformed to local coordinates, and 
the internal force vector is computed as described above. One thing to be noted in Eq. (35) 
is that if @ is initially 0, the yielding criteria siabie is identical to the von Mises yield | 
condition. 

After calculating the strain tensor using Eq. (23), any previous plastic strain is 


subtracted using 


Let fate }- {er} (36) 
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The stress is then calculated using Eq. (24) with the components of te : Using these 
values, the hydrostatic stress, deviatoric stress and effective stress are calculated as follows: 
p=%(o, +0, +0, (37) 


(s}={o} + p(s} 68) 


qz=f3(s? + 53 + 53 + 2(s3 + 53 +52) (39) 


where 6 is the Kronecker delta function: 

{5} = {11100}7 | (40) 
At this point, F is calculated from Eq. (35). If F is greater than zero, indicating plastic flow, 
iteration is required to determine the new porosity and change in plastic strain. The 
predictor-corrector method used in Aravas [2 1] is also used here. Using the values of p and 


g previously calculated as a first guess, correction factors are calculated using 








(cy $= [A]? {Al} ad 
where 
a ch 
Al 42 
Cs ee (42) 
Ke+e > -3GE+E 


[Al= 


(43) 








F OF FF 00 OF 0% OF u area 0p OF Oy 
+ As [KEE + gE 83 +E, aan; he, op TAF, 3G + ipo ate, [+ A&a Sedeg She, 


{cy I-14 | ay 
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and 











a 249 
ae 
ar 2 
"op 
” ee 
£=29,0 el sinh| 2? 
260 206 
2 
2 3q -3q,p 
gF=29,0 2 | cosh| —~2— 
ae a be 260 ) 
eof + 2g, DP sinh a8 
0 Oo Oo 
2 
; 3q —3q)P | 34) .. ,{ —39P 
oF_ =~ 2¢g@ 2. | cosh| —*2— | + —* sinh] —“= 
apea, ~~ EP 34, | Qo? 202 Qo? 
er _— 4 
Ogdo> o? 


The values for a and f are used to correct the change in strain caused by hydrostatic pressure 
and the change in strain caused by effective stress (See Ea. (46)), which are then used to 


determine the change in the plastic strain vector, as shown in Eq. (47). 


Agt= Ag! +a (Ae5=0) - 
Agi=Aet! + B (Ah=0) 
3 
(Ag? j=4Ac, {5} + sel ]f} (47) 


This change in plastic strain, {A e? is then added to the total plastic strain, and the new 
elastic strain is calculated using Eq. (36), and the stress vector is calculated again using Eq. 
(24). Next, the change in void content, or porosity, is calculated as: 


AD=A D growth + A Dnucleation (48) 


A® gown =(1-®) ef (49) 
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where @, is the volume fraction of void nucleating particles and ¢, and sy are the mean and 
standard deviation of a normal distribution of nucleation strain, as suggested by Chu and 


Needleman [15] and utilized by Aravas [21]. Then the effective plastic strain, and void 


content are updated with 


p 7 DAE, + GAEq 
Ag? = — 51 
e (1-® Joo ( ) 
elem me? + Ae? (52) 


where ¢? is the effective plastic strain from the previous time step. At this point, the change 
in yield stress due to strain hardening is calculated (see below). If either @ or fis greater 
than a predetermined tolerance, the process iterates beginning with Eq. (41). The tolerance 


used for all examples presented in the paper is 1.0 x 10%. _ 


B. IMPROVING STABILITY IN THE CONSTITUTIVE EQUATIONS 


Stability in the procedure outlined above-is very dependent on the order in which the 


"various equations are evaluated. Aravas discusses the stability problem when considering 


problems involving large plastic strains [21]. While the change in plastic strain from one 
time step to the next will theoretically be small, high strain rate and changes in the strain- 
hardening characteristics of the material act to degrade stability. Equations (42), (43) and 
(45) are not complete in that they do not contain all of the derivative terms required by the 


chain rule. 


20 











The predictor-corrector method used is based on Newton's method, where the 


correction factors, a and B, are found from 








Fac, “Sias, |{@ 
Ac, Ac, ” I (53) 
Src, Fre, B Fe 
where 
2 3 . 
f -(4 +24 0cosi{ LP }-(+40") (54) 
on 20 
fy=he, Hi +hé, Ty (55) 
” @q | * @p 
Pach Of, O® Of, 005 
= +b 
ius, = * a) * 30 dhe, * do, dbs, Ce 
Ticks = — Gh Hr, OP | Ai 00% (57) 
oy 0q OM O0AE, Oo, OAE, 
2 2 2 2 
fr, =Livns, OD — coh, Fh @ | PF a (58) 
= aa 00,09 OA, op" ” paw OAé, Opdo, OAe, 





f -a, Ae | Of 28 , PA oy |-36244 Of, do, (59) 
oes "| Qpe® BAe, Apo, Ade, Gq’ Agdo, dhe, 


Derivative terms that are zero have been dropped, and K and G are the bulk and shear 
moduli, respectively. f, is Gurson's void function, and f, 1s the flow rule, both of which are 
driven to zero. As void content increases, the number of iterations required to achieve © 
convergence increases dramatically, and often diverges, instead. By removing all partial 
derivatives relating to void content, the stability of the algorithm is greatly increased, at the 


cost of only 2 to 3 extra iterations, depending on the current void content. The dependence 
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on the current yield stress is retained to allow convergence across a transition in the slope of 
the yield stress versus strain relationship. 

The partial derivatives of yield stress with respect to 4g, and Aé, are computed by 
calculating the slope from the current yield stress to the yield stress corresponding to the 
increase in plastic strain from the corrected values of Ag, and Agé,. This means that these 
derivatives will be zero on the first iteration, since both of the adel variables are initialized 
to zero. Several error traps must be included to prevent round-off and truncation error from 
causing the algorithm to diverge, and to prevent porosity from decreasing or becoming 
negative. One of the assumptions used in this implementation is that once voids form, they 
do not disappear. In other words, voids do not disappear when the element is placed in 


compression. 


C. STRAIN HARDENING 


Modeling the nonlinear elasto-plastic behavior of the material used is simplified by 
constructing a piece-wise linear version of the stress-strain plot. Using the tangent modulus, 


E,, for each piece-wise region, the yield stress is calculated with 


| j-! 
of *=00 + S Ene - 61-1) + Eg lef - 2-1) (60) 


i=] 
where o/*”' is the yield stress for this time step, o> is the original yield stress, ¢ is the 
total effective strain for the time step under consideration, and ¢;is the upper strain limit of 


the ith linear segment. ¢, is the original yield strain, and is calculated by the program as 


Z2 








simply oo /E. The current effective elastic strain is calculated by dividing the current yield 
stress by the elastic modulus, E. This is added to the current effective plastic strain to obtain 
the total effective strain. 
As an example, let the total effective strain from Eq. (49) lie within the second work- 
hardening segment. The new yield stress is 
of =03+ En(e- 0) + En(efa - a1) (61) 
Figure 3 illustrates this example. The algorithm calculates the new yield stress as a function 
of the cumulative effective plastic strain, the current effective elastic strain, and the original 


yield stress for each iteration of the damage constitutive equations. 





eff ¢ 
Eo e Be. 2 € 


£3 


Figure 3. Calculating a New Yield Stress. 
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IV. HOURGLASS MODE CONTROL 





Only one integration point is used in each plane parallel to the reference plane, which 
results in under-integration in the ¢7 plane. This leads to spurious, hourglass, or zero- 
energy modes in the element, which will yield useless results if left uncontrolled. The effects 
of these modes are shown in Fig. 4, a pinched cylinder where one eighth of the structure is 
modeled by utilizing symmetry boundary conditions. Clearly, no useful information can be 


obtained from the resulting solution. 





Figure 4. Hourglass Modes in a Pinched Cylinder Model. 
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Belytschko, et al., proposed an efficient means of controlling the hourglass modes 
of a similar element [251 The method described uses a portion of a stiffness matrix 
generated by full integration in all directions to modify the stiffness matrix generated by 
under-integration. Although the formulation proposed in this paper does not use a stiffness 
matrix, a similar approach is just as effective in controlling these modes. 

Rather than fully integrating in all three directions, the element is fully integrated in 
the £7 plane, but under-integrated in the ¢ direction, as shown in Fig. 5. The procedure 
described above for calculating the internal force vector is followed to generate an internal 


force vector related to these four integration points. 





Figure 5. Hourglass Mode Control Integration Points. 


The algorithm for calculating strain, stress, and force for the hourglass mode control 


integration points is identical to the algorithm used to produce the main internal force vector, 
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' with the exception of plastic strain. The damage constitutive equations described in the 
previous section are not utilized for hourglass mode control. The internal forces generated 
_ from the two integration schemes are treated like the stiffness matrices in Belytschko et al. 


[25]. For nz integration points through the element thickness, the new force vector is 


{Fu S={ Fe" } + AL S sur (62) 
where 
Ae Crm ot ae ee (63) 
and 
h -fe (64) 


The variable h is used here in Eqs. (62) and (64) aia of the used in Belytschko 
et al. [25], to — confusion with the various strains discussed in this paper. The variables 
used to calculate / are the element thickness (f) and the element’s surface area (A). The 
effect of r follows that described in Belytschko et al. [25], and is set to 0.05. The range of 
values for 7 that effectively controls the hourglass modes, but does not greatly affect the 
overall element stiffness is roughly 0.046 to 0.057 (determined experimentally). Since the 
elements are in arbitrary orientation in 3-D space, the area calculation is computed as the 
sai of the area of the two triangles formed by dividing the element at the diagonal between 


nodes one and three: 


A=LJ/P2 - p? + BE - D? (65) 


where 
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Dy = (y ~ Xp ) Oey ~ Xf) + (B63 - Xp ) (2 -¥Q) + (5-3) 3 - 5) 
(66) 


Pa. whee 2S { c44eS> 2A 3 
Do= (%, - X ) Oy - X, ) + (X2 - x2) (0% - x) + (x3 - x3) (x; -x}) | 


and 





(67) 








and (xc, x, xt ) is the location of node & in the global coordinate system. For these 


calculations, the element is assumed to be flat (no curvature along either the ¢ or 7 — 
directions). The effectiveness of this method of control is shown in Fig. 6, the same problem 


as shown in Fig. 4, but with the hourglass mode control described above. 


Z axis 





y axis x axis 


Figure 6. Pinched Cylinder Model with Hourglass Mode Control. 
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All verification problems analyzed in the following section were completed with 
hourglass control enabled. The current implementation uses hourglass control consistently, 
but allows easy modification to make hourglass control a user-defined option. The internal 
hourglass mode control can be disabled when the element is used in a general finite element 


program that includes various methods of spurious mode control. 
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V. IMPLEMENTATION IN RESEARCH-ORIENTED FE CODE 


An in-house general purpose finite element analysis program, referred to herein as 
FEA, was used to verify the algorithms discussed. Several of these verification problems are 
presented in the next section. The program is written in FORTRAN, and is tailored for easy 
modification, rather than si or ability to handle very large problems. The program reads 
the system parameters from a single input file, then writes the results to two ASCII text files: 
one contains the nodal coordinates and displacements for all time steps, the other contains 
the data associated with each element for all time steps (stress and strain tensors, effective 


plastic strain, void content, yield stress, etc.). 


A. PREPROCESSOR 


Since this program. is locally generated, there is no preprocessor or postprocessor 
available. In order to facilitate running several problems, a preprocessor was written in 
MATLAB that generates the proper input file for seta simple geometric shapes. The 
MATLAB script file will create a properly meshed flat plate, open box, cylinder, complete 
spherical cap, or a spherical cap with a hole in its top. The listing for the preprocessor is 


presented in Appendix A. 


B. | POSTPROCESSOR 


A postprocessor, also written as a MATLAB script file, utilizes a graphical user 


interface (GUI) to create the required graphs. It also has the capability to create an animated 
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display of the displacement over time. The displacements are scaled to make them clearly 
visible. This feature is useful in determining the effectiveness of hourglass control and 
verifying the boundary conditions applied. However, the MATLAB movie is memory 
intensive and slow, and is therefore useful only as a qualitative verification of proper 
problem definition. All graphs of problem solutions, except as noted, were generated using 


this postprocessor. The full listing of all associated files is provided in Appendix B. 


C. FEA SHELL ELEMENT SUBROUTINE DESCRIPTION 


The FORTRAN source code for the FEA implementation of the element presented 
here, along with the proposed hourglass control and constitutive model, is listed without 
comment in Appendix C. Since the primary purpose of the FEA program is to assist in the 
development of finite element algorithms, this subroutine is written with readability, rather 
than computational ericiency, as a primary concern,. A brief outline of the subroutine is 
shown in Fig. 7. The FEA program treats each type of element (shell, beam, etc.) in a single 
block, and loops through each element, then each integration point within the element. 

The order of calculation when evaluating the constitutive equations is critical, and 
any variation will cause instability. The two control variables, dep and deq, are updated by 
the correction factors on each iteration. All other variables are initialized to the values from 
the previous time step within each iteration. This essential feature allows stability of the 


algorithm. As discussed above, the variation of void content with dep and deg is not used, 
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but the variation of yield stress with the two control variables is calculated and included in 


the full expansion of the partial derivatives. 


33 





Re ON Pe 


te 


Program Outline for ELSH4L.F 


Enter subroutine 

Set constants 

Load gauss points and weights 
Begin loop over each shell element 


A. Initialize element internal force vectors 
B. Retrieve element material properties 
C. Compute element indices into system vectors 
D. Load nodal coordinates and displacements 
E. Compute local coordinate system and surface area 
F. Compute rotation matrices 
G. Transform displacements to local coordinate system 
H. Calculate hourglass force vector 
I. Begin Integration Loop 
1) Load shape functions and derivatives - . 
2) Calculate Jacobian, its inverse, and its determinate 
3) Calculate strain-displacement matrix ([B]) 
4) Compute strain ({e} = [B] {d}) 
5) Transform strain to local coordinate system 
6) Subtract any previous plastic strain from strain tensor 
7) Load last values of 6, 6, and €? 
8) Compute stress ({o} = [D]{¢}) 
9) Compute yield function (F) 
10) If F > 0, set dep = deq = 0 and begin iteration 
a) Calculate required partial derivatives 
b) Prevent Numerical Errors (Divide by Zero, etc.) 
c) Retrieve Previous 9, €?.., and 0, 
d) Update dep and deq 
e) Compute Incremental Plastic Strain Tensor 
f) Compute Void Growth 
g) Compute New o,, $;, and Peg 
h) . Compute Void Nucleation Factor 
1) Compute New @ and op 
) Recalculate Yield Function (F) 
h) If F is not Zero, go to step a) 
11) Last values used are retained. 
12) Transform stress to global coordinate system 
13) Calculate internal force for this integration point and sum 
14) Store material tensors and constants in system vectors 
15) End of integration loop 
J. Apply hourglass correction force 
K. Transform internal force vector to global coordinate system 
L. Compute element mass vector 
M. Load mass and force vectors into system vector 
N. End of element loop 


End of subroutine 


Figure 7. Outline of FEA Implementation. 
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VI. NUMERICAL EXAMPLES 


The element presented here has been extensively tested using a variety of problems 
in which an analytic solution was available, and has produced satisfactory results in all cases 
studied. The transformation matrices and constitutive equations were verified with 
specifically tailored problems, and the element passes the patch test. The examples presented 
below are representative of the test cases used:to verify the element, and are presented in the 
order of curvature: flat, singly curved, and doubly curved. For each case, an elastic solution 


is compared to available results, and then the elasto-plastic solutions are presented. 


A. ELASTIC PLATE 


A plate clamped on all four sides 1s subjected to a concentrated force at its center. 
The elastic modulus is 10 Msi (68.95 GPa), the density is 0.1647 slugs/in’ (0.147 kg/cm’), 
and Poisson’s ratio is 0.2. The dimensions of the plate are 10 in x 10 in x 0.1 in thick 
(25.4cm x 25.4 cm x 0.254 cm). The yield stress is set high enough to ensure a completely 
elastic response. Two integration points through the thickness are used. The applied force’ 
is 40 lbf (177.9 N). The finite element mesh uses symmetry to model one quarter of the 
plate, with appropriate symmetry eee conditions applied. Both a 2x2 (4 element) and 
4x4 (16 element) mesh are used in the finite element analysis. The results for both meshes 


and the analytic solution are shown in Table 1. Using a four-element mesh, the new shell 
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element obtained a displacement within 3.27% of the analytic solution, and 0.82% using a 


16 element mesh. 


Table 1. Comparison of Results for Elastic Clamped Plate. 


Center Node Peak Displacement (i) 
“4.74x10° 
494x107 
490x107 














B. THICK CLAMPED PLATE UNDER PRESSURE LOAD IN ELASTO- 
PLASTIC REGION 


A thick steel plate, clamped on all four sides, is subjected to dynamic pressure load. 
The plate is 6m by 6m by 0.6m thick (for a 10:1 ratio). Table 2 shows the material 


properties for the structure. Table 3 shows the properties for the void model. 


Table 2. Material Properties of Clamped Plate. 














‘Blastic Modutus(@)—Si SCO COCSC~C“—*~‘~sSC“‘CNSNCSC#POC*d 
“Tangent Modulus (E) +(| STO OSC~—SSC“‘“‘CS)N™COCPA CS” 
Yield Stress Sw) C~iCIOSSC“‘“‘NOCOC#C‘(#’SWPATSCS 
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Table 3. Void Characteristics of Clamped Plate. 


[Model Constmtq ——SSSSSsS—<“C~C‘“s*S*S*~dSSC“‘(SSNNCOC“‘LSSSC“NSC“SSCSC“‘C‘“ 
"Model Constant 7 —=S~—~—“—~S—“—C—s—“‘“‘“‘<‘ YS 









One quarter of the plate is modeled with a three by three element mesh, for a total of 
nine elements, with appropriate symmetry boundary conditions applied, and is shown below 


in Fig. 8. 





. Figure 8. Nine Element Clamped Plate Mesh. 
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Four integration points are used through the thickness.of each element. The 
calculation time step is 10° seconds, while the data is plotted at 2x10“ second increments. 
The analysis terminates at 0.04 seconds. The plate is subjected to a uniformly distributed — 
pressure acting downwards, and includes the appropriate drilling moment. The pressure 
increases linearly from 0.0 Pa at the start to 80 MPa at 0.01 seconds, then remains constant 
for the duration of the analysis. Three cases were analyzed: no void effects, void growth 
effects only, and void growth and nucleation. Figure 9 illustrates the effect of voids on the 
effective stress versus the effective strain relationship in the top of one of the border 
elements. When void effects are not included, the Von-Mises Equivalent (VME) stress 
follows the yield stress in the plastic region (See Fig. 9a). 

Including void effects causes the yielding before the VME stress reaches the yield 
stress. As the void content increases with continued plastic flow, the VME stress falls farther 
below the yield stress, as shown in Fig. 9b (see Eq. (35)). Table 4 summarizes the results 
of the three cases. The most significant difference is in the transverse normal stress, 6... 
This difference resulted in a 2.0% increase in the effective plastic strain, and a 2.6% increase 
in the peak deflection of the center of the plate. | 

Another interesting point is that the transverse normal stress was compressive and 
constant throughout the thickness up until plastic flow, as assumed in the formulation, then 
varied as the stress in the bottom fiber decreased and became tensile. Figure 10 shows the 
variation of transverse normal stress through the thickness of the element at the time of peak 


stress. Eventually, the transverse normal stress varied from compressive at the top, where 
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the pressure was applied, and tensile at the bottom. The examples that follow will not 
include a separate analysis for void growth effects only, since the difference of including 
nucleation effects at the resulting small plastic strains obtained do not cause significantly 


different results. 
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Figure 9. Void Effects in a Clamped Plate with Pressure Loading: Top (a) - No Void 
Effects, Bottom (b) - Void Effects. 
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Figure 10. Transverse Normal Stress Variation through Shell Thickness: Clamped Plate with 


Pressure Loading and Void Effects. 
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Table 4. Summary of Results for Clamped Plate with Pressure Load. 


Peak Values for Element #3 No Void Effects Void Growth Growth and Nucleation 
Oy (Gpa) 1.4498 1.3789 1.3784 









00588 0.0597 0.0597 


6,, (GPa) 0.6561 0.6030 0.6073 


(Max Tension) 
ia Co 
(Max Compression) 


0.0540 0.0551 0.0551 
® (Porosity) 0.0355 0.0357 


Deflection (m) -0.3801 | -0.3897 -0.3898 
(Center Node) 


o,, (GPa) 1.4501 1.3603 1.3598 






C. THICK CLAMPED PLATE WITH CENTRAL POINT LOAD IN THE 
ELASTO-PLASTIC REGION 


The geometry and material properties of this example are identical to those used in 
the previous example. The pressure load has been replaced with a point load at the center 
node. Four a are studied: without void effects or a drilling moment, with void effects, 
with a drilling moment applied, and with both void effects and a drilling moment applied. 

All four sides are clamped, and the center node displacement is restricted in the x and y 
directions. Then, a DM equal to the applied force times ee, the plate thickness is 
applied for both the no-void and void cases. The time history of the stress components in the 
center element is shown in Fig. 11, and the relationship between porosity and effective 
plastic strain is illustrated in Fig. 12. The results for all six cases are summarized below in 


Table 5. 
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Figure 11. Stress Component Time History in Bottom Fiber: Clamped Plate with 
Concentrated Load, Void Effects, and Drilling Moment Applied. 
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Figure 12. Porosity versus Effective Plastic Strain in Bottom Fiber: Clamped Plate 
with Concentrated Load, Void Effects, and Drilling Moment Applied. 
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Table 5. Summary of Results for Clamped Plate Subjected to Point Force. 
Element No DM No DM Drilling Moment Drilling Moment 
(Max Tension) 
(Max Compression) 


0.0838 0.0880 0.0812 0.0840 
D Porosity) 0.0663 0.0601 


Deflection (m) -0.4413 -0.4544 -0.4350 -0.4451 
(Center Node) 


The drilling moment increases the effective stress on the fiber in compression, and 












decreases the effective stress on the fiber in iil The reduction in the tensile stress 
results in a reduction in the effective plastic strain. Including void effects decreases the 
VME stress on the fiber in tension, and slightly increases the VME stress on the fiber in 
compression. The effective plastic strain is also inital All of these cau indicate that 
this formulation correctly predicts, in a qualitative sense, the effects of drilling moments and 
void growth and nucleation. Including only void effects increased the effective plastic strain 
by 5.0%, while including only the drilling moment decreased the effective plastic strain by 
3.1%. Including both void effects and the drilling moment increased the effective plastic 


strain by 0.2%, indicating the importance of applying both effects together. 
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| O,, (GPa) -0.0035 -0.1802 -0.0035 -0.1722 
(Max Compression) 


D. SIMPLY SUPPORTED PLATE WITH CENTRAL POINT LOAD IN THE 
ELASTO-PLASTIC REGION 


The material properties, geometry, mesh, and time values from the clamped plate 
problem above are used here. The magnitude of the force is 8x10° N, and the drilling 
moment, when applied, is —2.4x10°N. The same four cases are analyzed: no void or 
drilling moment effects, drilling moment effects only, void effects only, and both void 


and drilling moment effects. The analysis results for all four cases are summarized in 


Table 6. 


Table 6. Summary of Results for a Simply Supported Plate with a Point Force. 


Element No DM Drilling Moment No DM Drilling Moment 


o,, (GPa) 0.0060 © =0.0003 0.0114 -0.0003 
(Max Tension) 


























gPlastic (effective) 0.1277 0.1242 0.1382 0.1327 ° 
® Porosity) 0.1197 0.1092 


Deflection (m) -0.9219 -0.8990 -0.9770 -0.9420 
(Center Node) 


The qualitative results for this example are the same as for the clamped plate. Since 





the applied force causes greater initial yielding, void growth and nucleation has a greater 
effect. The drilling moment, when added to the void effects, reduces the amount of effective 


plastic strain in tension and displacement. 
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E. ELASTIC PINCHED CYLINDER 


An open-ended cylinder of radius 5.0 in. (12.7 cm), length 10.35 in. (26.289 cm), and 
thickness 0.094 in. (0.23876 cm) is subjected to a pinching load of 100 Ibf (444.82 N) (See 
Fig. 6). The elastic modulus is 10.5 msi (72.4 Gpa), Poisson’s ratio is 0.3125, and the 
density is 3.125x10° slugs/in’ (2783 kg/m’). The load is applied as a step function beginning 
at time t= 0 seconds. Using symmetry, the problem was reduced to a one-eighth section of 
the cylinder. The dynamic value should be twice the analytic static value. Inextensional 
shell theory gives a static radial contraction of 0.1117 in. (0.28372 cm). The maximum 
radial contraction of the model is 0.1995 in. (0.5067 cm), which translates to a static radial 
contraction of 0.09975 in (0.2534 cm). Using a j56-eleceit mesh (16 by 16), the maximum 
radial contraction was 0.2207 in. (0.5606 cm), for a static contraction of 0.1104 in (0.2804 
cm). The results are summarized in Table 7. The 16-element solution is within 10.7% of 


the analytic solution, while the 256-element mesh is within 1.21%. 


Table 7. Comparison of Results for Elastic Pinched Cylinder. 
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F. THICK PINCHED CYLINDER IN THE ELASTO-PLASTIC REGION 


The eer and void properties shown in Tables 3 and 4 are used for an open-ended 
cylinder with a pinching force at its center. The top half of the cylinder is modeled using a 
36 element mesh with appropriate symmetry boundary conditions along the bottom edge of 
mesh, as shown in Fig. 13. The analysis is calculated in 10° second steps, with output every 
2x10 seconds, from 0.0 seconds to 0.04 seconds. The pinching force is 66.792 kN a each 
side, with a drilling moment of 848.258 Nm applied on the on indicated. The same four 
cases are compared: no voids or drilling moment, void effects only, drilling moment only, 
and both void effects and drilling moment (DM). The effects of void growth and nucleation 
and the drilling moment is shown in the stress component time histories; Fig. 14 shows the 
stresses eithoit voids or a drilling moment, and Fig. 15 shows the stresses with both effects 


included. The results for all cases are summarized in Table 8. 
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Table 8. Summary of Results for Elasto-Plastic Pinched Cylinder. 


Peak Values for Center No Voids No Voids Voids Voids and 
Element No DM Drilling Moment No DM Drilling Moment 
Con: /Gy,(MPa) i ssid | 262.40 | 40 | 265.05 05 | 259.28 28 | 283.49 49 


Eesrective (X 10°) l ar ae 5722 ara 4273 ae 7049 
182.46 160.38 170.12 150.31 


o,, (MPa) 291.71 283.77 291.83 297.97 
o,, (MPa) 5.0619 -0.0238 1.0270 7.5235 
(Max Tension) 


G,, (MPa) -7.6358 -28.739 -7.0364 -35.373 
(Max Compression) 


ePlaste (effective) (x 10°) 0. = 0.7583 0.5771 1.6923 
© Porosity) 10) [NA | _NA__-| 04773 [1970 


Deflection (mm) a 12.919 -0.0220 -0.0214 
Global Maximum 

Deflection (mm) ~31.750 -31.750 -31.750 -31.750 
(Center Node) alt | 

Global Minimum 


Since there is only a small amount of plastic strain, the effects of void growth and 













“ 














nucleation are greatly reduced. Even with this small amount of plastic flow, the effective 
plastic strain was increased by over 116 percent with the inclusion of both drilling moment 
and void effects, while the peak center node displacement was only increased by 


approximately 4.5 percent. 
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Figure 14. Stress Component Time History in Bottom Fiber: Pinched Cylinder, 


No Void Effects or Drilling Moment. 7 
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Figure 15. Stress Component Time History in Bottom Fiber: Pinched Cylinder, Void 
Effects and Drilling Moment Applied. 


G. ELASTIC SPHERICAL CAP WITH A CENTER HOLE 


The results of analysis with the shell element developed here are compared to results 
of the problem proposed in MacNeal and Harder [26]. The structure is a hemisphere of 


radius 10 units with a thickness of 0.04 units, and has an 18° hole cut in the center. Taking 
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advantage of the symmetry in the structure, only 4 of the structure is modeled. An eight-by- 
eight mesh is used, for a total of 64 elements, with four integration points through the 
thickness of each element. The mesh used is shown in Fig. 16. The material properties of 
the structure are: E = 6.825x10’, p = 0.001, and v = 0.3. No void or drilling moment effects 
are employed in this example. Opposing forces of magnitude 2.0 are applied at each 
quadrant: the force at node 73 is of magnitude —1.0 and parallel to the y-axis, and the force 
at node one is of magnitude 1.0 and parallel to the x-axis. To obtain a representative static 
response using a dynamic model, the applied force begins with 0.0 magnitude at time 0, and 
increases linearly with a rise time of 0.16 seconds to the specified value. A calculation time 
step of 1x10* seconds is used, with termination at 0.22 seconds. The maximum deflection 
of node one was 0.0929, where the theoretical value is 0.0940, for a normalized displacement 
of 0.988. This is within the range of the results listed in MacNeal and Harder from the 
QUAD2 and QUAD4 elements, which are considered to be accurate for this problem [26]. 


These results are summarized in Table 9. 


Table 9. Comparison of Results for Spherical Cap with Elastic Loading, and Results from 


MacNeal and Harder [26]. 
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Figure 16. Mesh Structure for Spherical Cap. 


H. SPHERICAL CAP WITH A CENTER HOLE, IMPACT LOADING 


Key and Hoff [11] used the previous problem to test their element with an impact 
(step) load. The structure and matrial properties are the same, and the mesh is the same as 
shown in Fig. 16. The load is applied at its maximum at t=0.0" seconds, rather than ramped 
up. Figure 17, which plots the displacement of node one obtained from both the shell 


element presented here and the element developed by Key and Hoff [11], shows that the 


wees meee adem TNA a we etme = = 


results of this element are comparable to those obtained by other elements under impact 


loading. 


Node #1 


Present sseersesecersee KOY & HOff _ 


Mean Node Displacement 





0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 
Time (sec) 


Figure 17. Node One Displacement Time History: Pinched Spherical Cap with Impact 
Loading. 


I. SPHERICAL CAP WITH A CENTER HOLE, ELASTO-PLASTIC 
LOADING 


The spherical cap structure shown in Fig. 16 is now used to verify the stability and 


convergence of the damage-constitutive equations under double curvature, and to illustrate 
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both the effects of void growth and nucleation and a drilling moment on a more complex 
structure. The material properties used are the same as shown in Tables 3 and 4. The 
structure has a radius of 5.0 meters and a thickness of 25 cm, for a radius to thickness ratio 
of 20:1, commonly considered the limit of thin-shell theory. Although a force is applied at 
each quadrant, all four loads are directed inwards. Therefore, only one load is applied at 
node 37 of the mesh shown in Fig. 16. The applied load is 5.9397x10° N, with a drilling 
moment of -7.4246x10° Nm applied for the cases indicated. A calculation time step of 
1x10° seconds is used, with output every 5x10“ seconds and stopping at 0.2 seconds. The 
load is initial zero, and increases linearly to its maximum at 0.01 seconds. The results shown 
in Table 10 are for the element nearest the applied load, where stress is maximum. The 
porosity versus effective plastic strain relationship is shown in Fig. 18. Due to the small 
amount of plastic strain, porosity is restricted to the linear zone, as shown. Figures 19a and 
19b illustrate the effective stress versus effective strain in the same element on the 
compressive side and tensile side, respectively. Note that as the structure returns from a peak 
displacement, the stress follows the elastic modulus, not the tangent modulus. This reflects 
the correct behavior of material: the elastic modulus is not affected by plastic deformation, 
only the yield stress is affected. The node displacement where the force is applied is shown 


in Fig. 20. This figure shows that there are at least two vibration modes in operation. 


> 





Table 10. Summary of Results for Elasto-Plastic Spherical Cap. 


Peak Values for No Voids No Voids Voids Voids and 



















S257 


Deflection (m) 0.0233 0.0246 0.0233 0.0246 
(Point of Loading) 








0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 
. Effective Piastic Strain . 10° 


Figure 18. Porosity versus Effective Plastic Strain in Inner Fiber: Spherical Cap with 
Void and Drilling Moment Effects. 
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Figure 19. Comparison of Void and Drilling Moment Effects in a) 
Compression (Top) and b) Tension (Bottom) for Spherical Cap. 
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Figure 20. Contact Node Mean Displacement Time History: Spherical Cap with Void and Drilling 
Moment Effects. | | 


During testing, it became apparent that this problem is not suitable for testing elasto- 
plastic analysis. Since the structure is concave, the entire structure quickly collapses once 
yielding begins. In addition, the single point used to prevent rigid-body motion also 
provided a stress concentration and additional yielding (the "corner" began folding over). 


This necessitated using a force that just causes plastic flow, but does not collapse the 
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structure or cause yielding near the anchored node. This is illustrated by an analysis of a5 
percent greater load than used for the analysis shown in Fig. 20. The resulting node 37 
displacement is shown in Fig. 21, and the deformed structure in Fig. 22. However, it is still 
apparent that the qualitative results obtained in the previous elasto-plastic examples carried 
through to this problem; both voids and the drilling moment decreased stress in fibers under 
tension, and increased stress in fibers under compression. Due to the small amount of plastic 
strain, the increase in porosity was extremely small, which minimized the effects of the 


voids. 
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Figure 21. Contact Node Displacement: Spherical Cap with 5% Greater Load. 
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Figure 22. Deformed Structure at End of Analysis: Spherical Cap with 5% Greater Load. 


60 








Vil. CONCLUSIONS AND RECOMMENDATIONS 

A new shell formulation was developed for transient dynamic analysis which 
includes both void effects with plastic deformation and the transverse normal stress with 
_ drilling moment. The element is compatible with most shell elements which have three 
translation and three rotation degrees of freedom per node. 

A numerically stable scheme was developed for Gurson's nonlinear constitutive 
equation model. The model includes both void nucleation and void growth. Furthermore, 
hourglass control was implemented into the algorithm to avoid spurious modes caused by 
the under-integration scheme. 

The drilling moment was important for thick plates and shells, It induced large 
transverse normal stress and affected the plastic deformation. Similarly, the effect of voids 
on plastic deformation was not negligible. Numerical studies show that the effective plastic 
strain increased up " fifteen percent due to the combined effects of voids and the transverse 
normal stress. | 

The shell element presented here was tested on many cases, some of which had 
| reference solutions for comparison. The ie element gave accurate results to those 
problems. 

Although the element formulation presented here peformed well in numerical 
studies, the problems studied are not in the area where this element is most likely to be 
useful. The number of published problems available for comparison is limited, especially 
in the case of dynamic plastic deformation. The solutions provided by this element in the 
plastic region need to be verified by experimental data, and the constants of Gurson’s void 
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model need to be determined for a variety of materials and porosity conditions. In all, there 
are six material property constants that must be determined before an accurate solution can 
be obtained: 9, 43» 43>. fus Cx» aNd Sy. These are not readily available, but play a vital role in 
the constitutive equations of this formulation. 

Implementing this shell element into a general purpose finite element program will 
allow easier comparisons with other element formulations, including solid brick elements, 
and will facilitate the analysis of more complex problems. This work has already begun for 
a version of DYNAS3D, but is not yet complete. 

| The method of hourglass mode control used implemented a full version of the same 
formulation used for the internal force calculation. Efficiency may be increased by using a 
simpler formulation to calculate the internal hourglass forces. Along the same lines, no 
attempt has been made to improve the efficiency of the code presented here so that 


readability may be preserved. 
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APPENDIX A. PREPROCESSOR IN MATLAB FOR FEA. 


% feapre.m - Preprocessor for fea program 

% Program Limitations: 

% No input of Initial Conditions 

S Boundary Conditions must be applied manually to result file 
S Load must be applied manually to result file 

% Structure will be homogenous 

% Only 1/4 cylinder supported at this time 


% Get Input data 
type = input('Type of Structure (l=cylinder, 2=plate, 3=shoe box, 
4=capl, 5=cap2):'); 


% Structure Specific Input 


nelx = input('Number of Elements along x-axis:'); 

nelp = input('Number of Elements perpendicular to x-axis:'); 
xl = input ("Lower X-axis value:'); 
x2 = input('Upper X-axis value:'); 

R = input ('Radius:'); 

thetal = input('Starting Angle (degrees):'); 

thetaz2 input ('Ending Angle (degress):'); 


elseif type== 
= input ("Number of Elements along x-axis:'); 
nelp = input('Number of Elements along y-axis:'); 


xl = input('Lower X-axis value:'); 
x2 = input('Upper X-axis value:'); 
yl = input ("Lower Y-axis value:'); 
y2 = input ("Upper Y-axis value:'); 

elseif type== 


nelx = input ("Number of Elements along x-axis:'); 


nelp = input('Number of Elements Perpendicular to x-axis:'); 
xl = input ("Lower X-axis value:'); 
x2 = input ('Upper X-axis value:'); 
yl = input('Lower Y-axis value:'); 
y2 = input ('Upper Y-axis value:'); 
zl = input ("Lower Z-axis value:'); 


z2 = input('Upper Z-axis value:'); 


elseif type== 
R = input ('Radius:'); 
alpha = input('Phi (angle from Z-axis of inscribed circle):'); 
neap = input('"Size (1=25 elem 2=42 elem 3=63 elem):'); 
if ncap== 
narc=6; 
ntheta=11; 
elseif ncap==2 
narc=8; 
ntheta=13; 
elseif ncap==3 
narc=10; 
ntheta=15; 
else 
errordlg('Bad Type"); 
end 
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elseif type== 

R = input ('Radius:'); 

alphal = input('Alphal (elevation angle to bottom Laing: (1n 
degrees):' 


° 
f 


) 
alpha2 = input('Alpha2 (elevation angle to top ring (in degrees)):'); 
thetal = input('Thetal (start-angle from x-axis (in degrees)):"); 
theta2 = input('Theta2 (end angle from x-axis (in degrees)):'); 


nelp = input('Number of Elements along elevation angle:'); 
nelx = input('Number of Elements along theta angle:'); 
else | 
errordlg('Bad Type'); 
end 


 & Input Material and Time Data 


E = input('Elastic Modulus:'); 

rho = input('Density (mass / unit volume):'); 

pois = input('Poissons Ratio:'); 

t = input ('Thickness:'); 

Syp = input('Yield Stress:'); : 

ipt = input ('New(2) or Old(0) shell element (enter 2 or 0):"'); 
stime input ('Start Time:'); 

etime input ('Stop Time:'); 

delt = input('Calculation Step Time (delta):'); 

ptime = input('Print Step Time:'); 


2 If old shell element, adjust density to mass / unit area 
if ipt== 
rho = rho * t; 
ipoint = 0; 
else 
ipoint = input('Number of Integration Points in Z direction:'); 
end 


% *** Cylinder *** 

if type== 
% Calculate Node x,y,z coordinates 
thetal = thetal * pi / 180; : 
theta2 = theta2 * pi / 180; 


nndsx = nelxti; 

nndsr = nelptl; 

xstep = (x2 - xl) / nelx; 
rstep = (theta2 - thetal) / nelp; 
i= 1; 


ne = nelx * nelp; 
mnodes = nndsx * nndsr; 


% Initialize Vectors 
x = zeros(nnodes,1); 
Y = xX; 

Zz 

nodes = zeros(ne,4); 


for xt=xl:xstep:x2 
for r=thetal:rstep:theta2 


x(i) = xt; 
y(i) = -R * cos(r); 
z{i) = R * sin(zx); 
i=i+i1; 
end 
end 
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xa 
yd 
za 


end 


i oi i 
NK x 


™ OS 
ve) 
~*~ 


% k*k* Plate Bec 


if type== 
nndsx = nelx+1; 
nndsy = nelptl; 
xstep = (x2 - x1)/nelx; 
ystep = (y2 - yl)/nelp; 
1 = Le 


ne = nelx * nelp; 
nnodes = nndsx * nndsy; 


xX = zeros (nnodes,1); 
Y = XX, 
Z= xX 


nodes = zeros(ne,4); 


end 


xd 
yd 
za 


zeros (nnodes,1); 
xd; 
-~1 * ones (nnodes,1); 


for xt = xl:xstep:x2 


for yt=yl:ystep:y2 


x(i) = xt; 
y(i) = yt; 
i =i + 1; 
end 
end 


% *** Shoe Box *** 
if type== 
mndsx = nelx + 1; 


nndsp = nelp + i; 

xstep = (x2 - x1) /nelx; 
ystep = (y2 - yl)/(nelp 
zstep = (z2 - z1)/(nelp 
i= 1; 


ne = nelx * nelp ; 

nnodes = nndsx * nndsp; 
= zeros (nnodes,1); 

— x; 


? 
7 
7 
nodes = zeros (ne, 4); 


for xt = xl:xstep:x2 
% Front 


for zt = zl:zstep: 


<CL). = "2s 
y(i) = yl; 
yda(i) = 1 


end 


Z2 
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% Top 
for yt = yltystep:ystep: y2 
x(i) = xt; 
y(i) = yt; 
z{i) = 22; 
za(i) = -1; 
; 1 aed 
| end 
3 %$ Back 
for zt = z2-zstep:-zstep:zl 
x(1) = xt; 
y(i) = y2; 
ya(i) = -l; 
Z(i) = 2t; 
Die es Ma Oa 
end 


end 


Ce, ee ee 


% Spherical Cap #1 
if type==4 
& Generates FEA mesh of spherical cap using symmetry 
% Also handles nodal connectivity 


% Input Values 
theta = linspace(0,pi/2,ntheta) ; 
alpha = pi * alpha / 180; 

phi = linspace(0,alpha,narc) ; 


% Spherical Coordinates for bottom nodes 
for i=l:ntheta 


P(i,:) = [R theta(i) phi(narc) ]; 
bet (i) = 7; 

ber(i) = 7; 

end 


% Spherical Coordinates for next row of nodes 
mnodes = ntheta; 
for i= 1: (ntheta-1)/2 
P(itntheta,:) = [R theta(2 * i) (phi (narc)-(phi(narc) -... 
phi (narc-1))/2)]; 
end 


% Spherical Coordinates for the rest of the nodes 
nnodes = nnodes + (ntheta-1)/2; 
-nrow = ntheta - (ntheta-1)/2; 
for j=narc-1:-1:1 
Th = 0; 
dTh = pi / (2 * (nrow - 1)); 
bet (nnodes+1)=2; 
ber (nnodes+1)=6; 
for i=l:nrow 
P(itnnodes,:) = [R Th phi(j)]; 
Th = Th + dtTh; 
end 
nnodes = nnodes + nrow,; 
bct (nnodes) =1; 
ber (nnodes) =5; 
nrow = nrow - 1; 
end 
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% Change BC of top node 
bet (nnodes) =4; 
ber (nnodes) =7; 


% Convert to Cartesian Coordinates 
for i=l:nnodes 


x(i) = P{i,1) * cos(P(i,2)) * sin(P{ 
y(i) = P(i,1) * sin(P(i,2)) * sin(P( 
z(i) = P(i,1l) * cos(P(i,3)); 

end 


% Element Connections 
%$ Bottom row first 
mid = ntheta; 

top = 








1,3)); 
1,3)); 


mid + ntheta - (ntheta - 1)/2 -1; 
ne = 1; | ; 
for i=1: (ntheta-1)/2, 
cl = (i-1)*2 +1; 
nodes(ne,:) = [cl (cl+1) (midti) (topti)]; 
nodes(netl,:) = [(cl+1) (cl+2) (topt+l+i) (mid+i)]; 
ne = ne + 2; 
end 
% Rest of elements 
nrow = ntheta - (ntheta-1)/2 - 1; 
bot = mid; 
mid = top; 
top = top + nrow + 1; 
for j=l:nrow 
for i=l:nrow 
cl = mid + i; 
nodes(ne,:) = [cl (bott+i) (cl+1) (topti)]; 
ne=netl; 
end 
bot=midt+1; 
mid=top; 
top=toptnrow; 
nrow=nrow-1; 
end 
ne = ne - 1; 
end 
% Spherical Cap #2 (Has Hole it top) 
if type==5 
nndsx = nelx+l; 
nndsy = nelp+l; 
theta2 = theta2 * pi / 180; 
thetal = thetal * pi / 180; 
alpha2 = alpha2 * pi / 180; 
alphal = alphal * pi / 180; 
tstep = (theta2 - thetal) /nelx; 
astep = (alpha2 - alphal)/nelp; 
i= 1; 
ne = nelx * nelp; 
~nnodes = nndsx * nndsy; 
x = zeros(nnodes,1); 
yY = &; 
ZS 
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nodes = zeros(ne,4); 


for thetat = thetal:tstep:theta2 
for alphat=alphal:astep:alphaz2 
x(i) = R * cos(alphat) * cos(thetat); 
= R * cos(alphat) * sin(thetat); 
z(i) = R * sin(alphat); 
i= i #1; 


end 


end 


% Calculate Element Connectivity (will work for all shapes except cap) 


if 


type~=4 
for i=0:(nelx - 1) 
for j=l:nelp 
ie = (i * nelp) + j; 
nodes (ie,1) (i * (nelp + 1)) + 3; 
nodes (ie, 2) nodes(ie,1) + nelp + 1; 
nodes (ie, 3) nodes(ie,2) + 1; 
nodes (ie, 4) nodes(ie,1) + 1; 


end 
end 


end 


€ Void content 
PhiO = input('Initial Void Content:'); 


input ('Value for qli:'); - 
input ('Value for q2 (Enter 0 for Default):'); 

input ('Value for q3 (Enter 0 for Default):'); 

input ('Nucleating Particle Content:'); 

input ('Mean Nucleation Strain (Enter 0 for Default):'); 

input ("Nucleation Standard Deviation (Enter 0 for Default):'); 


loud ft a tt 


% Create Stress-Strain curve 
correct = 'N'; 
while ((correct == 'n') | (correct == 'N')) 


iseg = input('Number of Plastic Segments (0-9):"); 
if (iseg ~= 0) 
for i=l:iseg 
slope (i)=input ('Slope of Section:'); 
strain(i)=input('Right Strain limit of Section:'); 
end 
end 


strn(1) = 0; 
strss{1) = 0; 
strn(2) = Syp/E; 
strss(2) = Syp; 
if (iseg~=0) 
for i=l:iseg 
strn(2+i)=strain(i); 
strss(2+i)=strss(2+i-1) + (strn(2+i)-strn(2+i-1))*slope (i); 
end 
end 


figure (1); 
plot (strn,strss),grid,xlabel('Strain'),ylabel('Stress'),title('Mat 


erial Yield Property') 


correct = input('Is the Stress-Strain Plot Correct? (y or Dye Ss") 


68 








end 


$ Fill eprop 
eprop=zeros (1,40); 


eprop(1) = rho; 
eprop(2) = E; 
eprop(3) = pois; 
eprop(4) = t; 
eprop(5) = Syp; 
eprop(6) = ql; 
eprop(7) = q2; 
eprop(8) = q3; 
eprop(9) = £n; 
eprop(10) = en; 
eprop(1ll) = sn; 
eprop(12) = Phid; 
eprop(13) = iseg; 


if (iseg > Q) 
for i=l:iseg 
eprop (12+2*i)=slope (i); 
eprop (12+2*i+1)=strain (i); 
end 
end 


%$ Display resulting structure 
figure (2) | 
AdjView = [-37.5 30]; 
shownode; 

rotate3d on; 


% Apply Boundary Conditions 


if type~=4 
bct = zeros(l1,nnodes) ; 
ber = bct; 

end | 


disp('Boundary Conditions:'); 


disp('0 = No Constraint'); 
disp('l = Constrained in X'); 
disp('2 = Constrained in Y"'); 
disp('3 = Constrained in 2Z'); 
disp('4 = Constrained in X and Y'); 
disp('5 = Constrained in Y and Z'); 
disp('6 = Constrained in X and Z'); 
disp('7 = Constrained in X,Y, and Z"'); 
C180 aka) > 

if type==4, 


disp('({It is inadvisable to change BC''s for Spherical Cap)'); 
end 


nod = input('Node to change Boundary Conditions (0 to end, -1 to 
1ist}:*)3 
while ((nod ~= 0) & (nod <= nnodes) ) 
if nod < 0 : 
disp ('Node Trans Rot'); 
for i=l:nnodes 


Gisp(sprintf(' %d td $a',i,bet(i),ber(i))); 
end 
else ; 
disp(sprintf('Current: Translation = %d, Rotation = ',... 
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'S$da', bct (nod) ,bcr (nod) )); 


bet (nod) = input('New Translation Boundary Condition(0-7):'); 
ber (nod) = input('New Rotational Boundary Condition(0-7):"); 
disp(* '); 

end 


nod = input ('Node to change Boundary Condition (0 to end):'"); 
end 


% Create Load History 
disp (" (Pressure Load only works for Flat Plate)'); 


nload = input('Number of nodes to apply load to (0 = pressure load): 


ntime = input('Number of time intervals for Load Description:"); 
disp('Force Direction:'); 


disp('l = Force along X-axis'); 
disp('2 = Force along Y-axis'); 
disp('3 = Force along Z-axis'); 
disp('4 = Moment about X-axis'); 
disp('5 = Moment about Y-axis'); 
disp('6 = Moment about Z-axis'); 
if nload == 


nload = nnodes * 6; 
disp('Input Pressure Time History'); 
for i=l:ntime 


ldtime(i) = input('Time:'); 
Pressure(i) = input (° Magnitude(use neg. values for Beebo, 
"pressure):'); 
end 


F = zeros (nnodes,1); 
for elem=l1:ne 
% Get Coordinates of Nodes for this syement 


xl = x(nodes(elem,1)); 
x2 = x(nodes(elem,2Z)); 
x3 = x(nodes(elem, 3)); 
x4 = x(nodes(elem, 4)); 
yl = y(nodes(elem,1)); 
y2 = y(nodes (elem, 2)); 
y3 = y(nodes (elem, 3) ); 
y4 = y(nodes (elem, 4)); 
zi = z(nodes(elem,1)); 
z2 = z(nodes(elem,2)); 
z3 = z(nodes (elem, 3) ); 
z4 = z(nodes(elem,4)); 


% Calculate length of ides; then area 


Li = sqrt ((x2-x1) *2+ (y2-y1) *2+ (z2-z1)%2); 

L2 = sgrt( (x3-x2) *2+ (y3-y2) “2+ (z3-z2)%2); 

L3 = sqrt ((*x4-x3) *2+ (y4—-y3) *2+ (24-23) %2); 

L4 = sqrt ((x1-x4) *2+ (yl-y4) *2+ (z1-z4)%2); 

dotl = (x1i-x2)* (x3-x2)+(yl-y2) * (y3-y2) + (z1-z2) * (z3- Z2)}3 
dot2 = (xi-x4) * (x3-x4)+(yl-y4) * (y3-y4)+(z1- z4)* (z3-2z4); 
Area = 


0.5 * (sqrt(L1*2 * L2*2 - dot1%2) + sqrt (L3%2 
L4*%2 - dot2%2)); | 


ee 


% Calculate Force per unit pressure (assumes square elements) 


F(nodes(elem,1)) + Area / 4; 
F(nodes(elem,2)) + Area / 4; 
F(nodes(elem,3)) + Area / 4; 
F(nodes(elem,4)) + Area / 4; 


F (nodes (elem, 1) ) 

F(nodes (elem, 2) ). 

F (nodes (elem, 3) ) 

F (nodes (elem, 4) ) 
end 


nid ul 
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% Create Node, Force,Time Vectors 
ne = 1; 
for n=l:nnodes, 

loadnode(nc) = n; 


loadnode(nctl) = n; 
loadnode (nec+2) =n; 
loadnode (nc+3) = n; 
loadnode (nct4) = n; 
loadnode(nc+5) = n; 
loaddir(nc) = 1; 

loaddir(nct+l1) = 2; 
loaddir(nct2) = 3; 
loaddir(nc+3) = 4; 
loaddir(nct+4) = 5; 
loaddir(nect+5) = 6; 


loadtime (nc,:) = ldtime; 


loadtime (nct1,:) = ldtime; 
loadtime (nc+2,:) = ldtime; 
loadtime (nc+3,:) = ldtime; 
loadtime (nct4,:) = ldtime; 
loadtime (nc+5,:) = ldtime; 


xd(n) * F(n) * Pressure; 
yd(n) * F(n) * Pressure; 
zd(n) * F(n) * Pressure; 
loadmag(nc,:) * t / 2; 
loadmag(nctl,:) * t / 2; 
loadmag(nc+2,:) * t / 2; 


loadmag(nc,:) = 
loadmag (nct1,:) 
loadmag (nct+2, :) 
loadmag (nct3,:) 
loadmag (nct4,:) 
Lloadmag (nct5, :) 


nce = nc + 6; 


else | 
disp('Input Load/Time/Magnitude History'); 
for n=i:nload 

loadnode(n) = input('Node Number:'); 


loaddir(n) = input('Load Direction:"'); 
for nt = l:ntime , 
loadtime (n,nt) = input('Time:'); 
loadmag(n,nt) = input ('Magnitude:'); 
end a 
end 


end 


% Create draft input file 
F = fopen('input.in','w'); 


fprintf(F,'%d 1 $d td 0 $d O \r\n',nnodes,nload,ntime,ipt); 

fprintf(F,'%te te te te \r\n',stime, etime,delt,ptime) ; 

fprintf(F,'0 00000 \r\n'); 

fprintf(F,'0 $d 0 0 0 \r\n',ne); 

fprintf(F,'0 0 0 0 \r\n'); 

‘fprintf(F,'0 O\r\nO O\r\n%d\r\nl1 1 \r\n',ipoint); 

for 1=1:4 : 

fprintf(F,'te te te te Se te Se Fe Fe F$e\r\n',.. 

eprop ((1+(i-1)*10) :(i*10))); 

end 7 

for i=l:nnodes 

fprintf(F,'td td %d te %e $e \c\n',... 

i,bet(i),ber(i),x(i),y(i),z(i)); 
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end 


+ ae cert wees sete ete ee 


for i=l:ne 
fprintf(F,'td 1 %d $d $d $d O \r\n',i,nodes(i,:)); 
end 
for i=l:nload 
fprintf(F,'$d d\r\n',loadnode(i),loaddir(i)); 
for nt = l:ntime 
fprintf(F,'te e\r\n', loadtime (i,nt),loadmag(i,nt)); 
end 
end 


fclose(F); 
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APPENDIX B. POST-PROCESSOR IN MATLAB FOR FEA. 


% fea_disp.m 

%$ P. McDermott 

3 

% Uses GUI to read and display output from fea program . 
SESSEEEEEESESEEETEESSSSSSSESSSESEEEEESTEESSEESSSESS 


Sern = get(0, 'ScreenSize'); 
= round(0.9 * Sern(4) / 2 - 50); 
wd = round(0.9 * (Sern(3) - 120) / 2); 


figure ('NumberTitle','off','Name','FEA', 'MenuBar', 'none',... 
"Position', [20 (Scrn(4) - 500) 120 460], 'Color','b'); 


> 
Fh 
fl 


hfs = figure('NumberTitle', 'off', 'Name', ‘Structure 
View', 'Visib', 'off',... 
"Position', [160 (Scrn(4)-70-ht) wd ht]); 


hfd = figure('NumberTitle', 'off', 'Name', 'Displacement','Visib','off', 
"Position", [(wd+180) (Scrn(4) - 70-ht) wd ht]); 


hlast 0; 

AdjView = [-37.5 30]; 
flagd = 1; : 
flagv = 0; 

flaga = 0; 

flags = 0; 


File = ‘output.fea'; 
File2 = ‘output.fes'; 


P = pwd; 
Path = strcat(P,'\' v 


he change = uicontrol(hf, 'Style', 'push', 'String', 'Change File', 
'Position’, [20 380.80: 20] > ««-. 
“CaliBack’ ; [sas 
'[Filel Pathl]) = uigetfile(''*.fea'',''Load FEA Output 
|e ee ee ; 
‘if Filel ~= 0,'... 
"File = Filel;'... 
"Path = Pathi;'... 
'FileZ = File;:’ s.«-. 
'File2(length(File2)) = ''s'';'... 
"set (he_ file, ''String''.File);' 
Tend']); 


he_load = uicontrol (hf, 'Style', 'check', 'String','Load File',... 
"Position', [20 420 80 20],... 
"CallBack',[... 
"FilePath = [Path File];',... 
'FilePath2 = Each File2];',... 


"parse; ' 
‘set (hc_ clear, ''Enable'', lon! ye 
"set (hc disp,''Enable'',''on'');' 


‘set (hc_velo, ''Enable'',''on'');' 
"set (hc_acc,''Enable'',''on'');' , 
"set (nc_show, ''Enable'',"'on'');',... 
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"set (hc _load,''Value'',1,'"Enable'',''off'');']); 


he clear = uicontrol (hf, 'Style', 'push', 'String', ‘Clear 
Data', 'Enable','off',. 
'Position' [20 60.80: 20], * CallBack’, 1. 
‘clear n* x* y* z* 1* s* t* v* a* fp factor d etime',... 
COC: Me CIM ye 
‘set (hta,: *Virsib® \," "Ori" ) se" 724% 
"if (flagv>0),set(hfv,''Visib'',''off'');end;', 
'if (flaga>0),set(hfa,''Visib'',''off'');end;' 
set (hfs; 'Visib’';*  OLt! *)y* Fees 
'flagd = 1; flagv = 0; flaga = 0; flags = O;', 


"set (hc_load,''Value'',0, "'"Enable'', ao) 2 Maal ae 
‘set (hc _disp,''Enable'', WI OLE. * 5p Gece sd 

"set (hc_ ~velo,''Enable'', Toft" Ji pales 

'set (he acc,'’'Enable'',''off'');',... 

‘set (hc_show, ''Enable'’, Mofett);! 

"set (hc_replay,''Enable'', MTOEET");, 


‘hlast = 0;',... 
"set (hc_mview, ''Enable'', VOLPE eye x 
'set (hc clear,''Enable'', UT OEE ye dee 


he close = uicontrol (hf, 'Style', 'pusn ', ‘String! ») CLOSC 7 sad 

"Callback', [. 
(elose (hfdie 1.044 
"if (flaga>0),close(hfa);end;',... 
'af (flagv>0),close(hfv);end;',... 
"close (hfs);' 
“cléar;:’ 
ee eee 

'"Position', [20 20 80 20]); 


he_disp = uicontrol (hf, 'Style','push', 'String', 'Displace', 'Position"',... 
[20 260 80 20],... 
‘Enable', 'off', 
‘Callback’; [«.« 
"if (flags)', 
"'AdjView = get (hax, ''View'');', 
gi — 1 9 (ok Muara 
'figure(hfd);',.-. 
'displace;',... ! 
‘set (hc replay, ''Enable'',''on'');', 
'flagd 
"hlast 


ff ea 
hfd;']); 
he mview = uicontrol (hf, 'Style','push','String', ‘Adjustment’, 
"Position', {20 100 80 20],... 
‘Enable’, 'off', 
'CallBack',[... 
'figure(hfs);', 
'mmview3d;']); 


he file = uicontrol (hf, 'Style', 'text', 'Position', [20 340 80 20],... 
"String"; File); 
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he velo = uicontrol (hf, 'Style', 'push', 'String', 'Graphs','Position', [20 
220 80 20),... | 
"Enable', 'off',... 
"Callback',{[... 
TE LAOVE 1 2 yp ean 
*control;']); 


he_ace = uicontrol(hf, 'Style', 'push', 'String', 'Stress', 'Position', [20 
180 80 20],... 
“Pnaple: 7 OLt pues 
‘Callback’ , («2 
tat <(tlags) 75-5 
"AdjView = get (hax, ''View'');',... 
SONICS © 5.6 4-5 
"figure (nid) 7" 44- 
TStresscoLlLorys p44 
'tlaqd: = 07" ;.405 | 
"set (hc_replay, ''Enable'',*'on'');',... 
‘hlast = hfd;'j); 


he replay = uicontrol (hf, 'Style','push', 'String', 'Replay',... 
"Position', [20 140 80 20],... 
‘Enable’, 'off',... 
"Callback',[... 
"figure (hlast);',... 
"movie(m);']); 


he show = uicontrol(hf, 'Style', "push', 'String', 'Show',... 
"Position', [20 300 80 20],... © 
'Rnable', 'off',... 
"Callback',[.... 
"ticure (his) 7" 4s. 
"shownode;',... 
"set (hc_mview,''Enable'',"'on'');',... 
"flags = 1;',... 
"hax =.gca;']); 


SESESSESESESSESESEESTS End FEA DISP.M SS BS3SBsSsssssssssssSsSSSTBSBSFFs 
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% PARSE.M 

% This m-file reads the output file generated by fea.f 

% and parses it into matlab usable variables which will be used 
S by another program to graphically display the data in various 
g formats | 

% McDermott 1999 


SELYSSLLFSYSFSSALLLSELESLESESSESEEEEEETEEESEEEEESESSESESESSTTTESS 


% Create Progress Box and lst Message 
hprog = figure('NumberTitle','off', 'Name', 'Load File 
Progress', 'MenuBar', 'none',... 
‘Position’, [(Scrn(3)/2 - 80) (Scrn(4)/2 - 40) 260 80],'Color','y'); 
hmsg = uicontrol('Style', 'text', 'Position',... 
[20 45 220 20], 'BackgroundColor',"y',... 
'String', 'Opening Files...'); 
hmsg2 = uicontrol('Style', 'text', 'Position',... 
[20 15 220 20], 'BackgroundColor','y',... 
'String',FilePath); | 
drawnow; 


fp=fopen (FilePath, 'r"); 
fp2 = fopen(FilePath2,'r'); 


top=fscanf (fp, 'se', 33); 


% Update Progress Box 
set (hmsg, 'String', 'Loading Displacement"); 
drawnow; 


% split into var names 


set (hmsg2, 'String', "Input & eprop'); 
drawnow; 


nnodes=top (1); 
nmat=top (2); 
nload=top (3); 
loadtime=top (4); 
initc=top (5); 
stime=top (8); 
etime=top (9); 
dt=top(11); 
ne=top (19); 
ipoint = top(31l); 


for i=l:nmat, 
eprop = fscanf(fp, 'te',40); 
end; 


for i=1: (nnodes), 


fscanf (fp, 'te',3); % dummy index 
x(i)=fscanf (fp, '%e',1); % x coordinates of nodes 
y(i)=fscanf (fp, '%e',1); % y coordinates of nodes 
z(i)=fscanf(fp,'se',1); % z coordinates of nodes 
end; : 
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for i=l:ne, 
fscanf (fp, '%e',2); 


nodes (i,:) 


=fscanf (fp, '%e',4)'; 





fscanf(fp,'%e',1); 


end; 


initial conditions 
fscanf (fp, 'se', via (6eiedes) + (2 <nioad* Tisadeimasiy yy): 
nsteps=round ( (etime-stime) /dt) ; 


time L); 
temp nsteps; 
for t=l:nsteps, 


set (hmsg2, 'String', ['Time Step: 


drawnow; 


test=fscanf (fp, '%te',1); 
if (test > 0.0), 
time(t) = test; 
for i=l:nnodes, 


fscanf (fp, '%e',1); 


xu(i,t)=fscanf(fp,'%e',1); 
xv(i,t)=fscanf(fp,'%e',1); 
xa(i,t)=fscanf (fp, '%e',1); 
fscanf (fp, 'e',1); 

yu(i,t)=fscanf (fp, '%e',1); 


yv(i,t)=fscanf(fp,'%e',1); 
ya(i,t)=fscanf (fp, '%e',1); 


fscanf(fp, 'se',1); 


zu(i,t)=fscanf(fp,'%e',1); 


zv(i,t)=fscanf(fp,'%e',1); 
za(i,t)=fscanf (fp, '%e',1); 


fscanf (fp, '%e" ei 

thxu(i,t)=fscanf (fp, '%e',1); 
thxv(i,t)=fscanf(fp, 'te',1); 
thxa(i,t)=fscanf(fp, '%te',1); 
fscanf (fp, 'e',1); 

thyu(i,t)=fscanf (fp, 'te',1); 
thyv(i,t)=fscanf(fp, 'te',1); 
thya(i,t)=fscanf (fp, '%e',1); 
fscanf (fp, '%e',1); 

thzu(i,t)=fscanf (fp, '%e',1); 
thzv(i,t)=fscanf (fp, '%e',1); 


end; 
else 


thza(i,t)=fscanf (fp, 'te',1); 


" num2Zstr(t) 


% associate nodes w/ elements 


a ee num2str(nsteps) ]); 


dummy index 


% EOF reached, terminate loop storing maximum number of time steps 


t = nsteps; 
end; 
end; 


nsteps = size(zu,2); 


fclose(fp); 


% Update Progress 


box 


set (hmsg, 'String' ee eee Scales' ); 


set (hmsg2, 'String', 
drawnow; 


ue 


% Autoscale vectors 


ax 
dy 


abs (max(x) - 
abs (max(y) - 


i ot 


min(x)); 
min(y)); 


T7 


dz = abs(max(z) - min(z)); 


% Check for flat structures 
if dx == 


ax = 1; 
end 
if dy == 0 
dy = 1; 
end 
if dz == 0 
az = 1; 
end 
scale(1) = dx / max(max(abs (xu) 
scale(2) = dy / max(max (abs (yu) 
scale(3) = dz / max(max (abs (zu) 


dscale = min(scale) * 0.25; 
if dscale < l 
dscale = 1; 


end 

xu = xu * dscale; 

yu = yu * dscale; 

zu = zu * dscale; 

scale(1) = dx / max(max (abs (xv) 
‘scale(2) = dy / max(max (abs (yv) 
scale(3) = dz / max(max (abs (zv) 


vscale = min(scale) * 0.75; 
if vscale < l 

vscale = 1; 
end 


xv * vscale; 
yv * vscale; 
gv * vscale; 


. xv 


Vv 
ZV 


dx / max (max (abs (xa) 
dy / max(max (abs (ya) 
dz / max (max (abs (za) 


scale(l1) 
scale (2) 
scale (3) 


ascale = min(scale) * 0.5; 
if ascale < l 
ascale = 1; 


end 

xa = xa * ascale; | 
‘ya = ya * ascale; 
Za = za * ascale; 


)); 
Pe. 
))3 


e 
? 
* 
4 


f 


)) 
)) 
)) 


))3 
))3 
))? 


. 


78 











% Allocate arrays for stress-strain data 
stress = zeros (ipoint*ne*6,nsteps) ; 
strain = stress; 

plastrn = stress; 

void = zeros (ipoint*ne,nsteps) ; 
yieldstrs = void; 

effplstrn = void; 

vmstress = void; 

effstrain = void; 

estrain = void; 

$pstrain = void; 


%¢ Update Progress Box , 
set (hmsg, 'String', ‘Loading Stress & Strain'); 
drawnow; 


% Read Stress-Strain Data 

for t=l:nsteps, 
test = fscanf(fp2,'%s',1); 
set (hmsg2, 'String', ['Time Step: ' num2str(t) ' / ' num2str(nsteps)]); 
drawnow; 


if (test > 0.0), 
tl=fscanf(fp2, 'te',1); 
for elem=l1:ne, | 
fscanf (fp2, '%s',1); 
fscanf(fp2,'%d',1); 
for point=l:ipoint 
index = (elem-1)*6*ipoint + (point-—1)*6; 
fscanf (fp2,'%d',1); 
for 1. = 1<+6,; 
ind = i+ index; 
stress(i + index,t)=fscanf(fp2,'%te',1); 
strain(i + index,t)=fscanf(fp2, '%e',1); 
plastrn(i + index,t)=fscanf(fp2,'te',1); 
end 
indexl = (elem-1)*ipoint + point; 
yieldstrs(indexl,t) = fscanf(fp2,'%e',1); 
void(indexl,t) = fscanf(fp2,'%e',1); 
effplstrn(indexi,t) = fscanf(fp2,'%e',1); 


end 
end 
else 
t = nsteps; 
end 
end 


% Calculate Von-Misses Stress and an equivalent total strain value 


% Update Progress Box 

set (hmsg, 'String', ‘Calculating Effective'); 
set (hmsg2, 'String', ‘Stress and Strain...'); 
drawnow; 

fclose(fp2); 


% Adjust for total strain 
$totalstrn = strain + plastrn; 
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$ Get Effective stress and strain (doesn't work with multiple materials) 
for elem=l:ne, 
for point=1l:ipoint, 
index = (elem-1)*6*ipoint + (point-1)*6 + 1; 
index1 = (elem-1)*ipoint + point; 
vmstress (indexl1,:)=sqrt (0.5*( (stress (index,:) -. 
stress (indext1,:)).*2 + ... 
( stress(indext1,:) - stress(indext+2,:)).%2 t+... 
( stress(index+2,:) - stress(index,:)).*%2 + 6 *.. 
(stress (index+3,:).*2 + stress (indext4,:).%2 .. 
+ stress (indext5,:).%*2))); 


estrain(indexl, :)=sqrt (2* ((strain (index, :)-strain(index+1,:)).“%2.. 
+(strain(index+l1,:)-strain(indext+2)).*2 + (strain (index+2)... 
- strain(index,:)).*2) + 3*(strain(indext+3).*2 + .. 

strain (indext4).*2 + strain(indext5) .*2))/(2* (l+eprop(3))); 
end 
end 


close (hprog) ; 


SELTEFLESLSESEEESSFESB_~O-EMG of PARSE.M SSSESELESSBSBTSSTSVBFISEBBEBBSB 
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% Control.m | 

% Display Von-Misses Stress and Yield Stress curves 

% | McDermott 1999 
SSSESESEEESEEESETSESELESEEEEEESESESEEETEETESESESSESEEETEEESESESEEESESESSESES 


hfv = figure('NumberTitle', 'off',"Name', 'Graph Control','Visib','on',... 
"MenuBar', 'none', 'Color','b', 'Position', [160 50 wd ht]); 
hfa = figure('NumberTitle', 'off','Name','Plot Window', 'Visib','off',... 


"Position', [(wd+180) 50 wd ht)); 


elementstr = ['l|']; 
for i=2:ne 
elementstr 
elementstr 
end 


strcat (elementstr,num2str(i)); 
strcat(elementstr,'|'); 


pointstr = ['1l]']J; 

for i=2:ipoint 
pointstr = strcat(pointstr,num2str(i)); 
pointstr strcat (pointstr,'|"); 

end 


nodestr = ['l|']; 

for i=2:nnodes 
nodestr = strcat (nodestr,num2str(i)); 
nodestr = strceat(nodestr,'|'); 

end 


typestr = ['Stress vs Strain|Stress vs Time|Strain vs time]Porosity vs 
Time | 4s. | 
"Porosity vs Strain|Sigma zz vs Thickness|Sigma xx vs Time|',... 
"Sigma yy vs Time|Sigma zz vs Time|Sigma Multiplot|',... 
'X Disp. vs TimelY Disp. vs Time|Z Disp. vs Time|',... 
‘Theta x vs Time|Theta y vs Time|',... . 
‘Theta z vs Time|MeanDisp’ vs Time|Plastic Strain|Elastic Strain']; 


% Start at Bottom 

hv_plot = uicontrol (hfv, 'Style', 'push', 'Position',... 
[(wa/2 - 40) 20 80 20], 'String', 'Plot', 'Callback’,... 
'flaga=1;drawstress;'); | 


hv_typetxt = uicontrol (hfv, 'Style', 'text', 'Position',[ 20 60 120 20],... 
"String', 'Graph Type:'); : : 


hv_ type = uicontrol (hfv, 'Style', 'popupmenu', 'Position',... 
[ 160 60 160 20 j],'String',typestr, 'Value',1,... 
‘Callback’, 'type = get(hv_type, ''Value'');'); 
type = 1; 


hv_pointtxt = uicontrol (hfv, 'Style', 'text", 'Position',... 
[ 20 100 120 20 J], 'String', 'Integration Point:'); 


hv_ point = uicontrol (hfv, 'Style', 'popupmenu', 'Position',... 
[ 280 100 40 20 ],'String',pointstr, 'Value',1,... 
7 "Callback', ‘point = get(hv_point,''Value'');"); 
point = 1; 


hv_elemtxt = uicontrol(hfv, 'Style', 'text', 'Position"’,... 
[ 20 140 120 20 ],'String','Element Number:'); 


$1 


hv_elem = uicontrol (hfv, 'Style', 'popupmenu', 'Position',... 
[280 140 40 20 },'String',elementstr, 'Value',1,... 
'Callback', 'selem = get (hv_elem, ''Value'');"'); 
selem = 1; 


hv_nodetxt = uicontrol (hfv, 'Style', 'text', 'Position',.. 
[20 180 120 20 ],'String', 'Node Number:'); 


hv_node = uicontrol (hfv, 'Style', 'popupmenu', 'Position',... 
[280 180 40 20], 'String',nodestr, 'Value',1,... 
'Callback', 'snode = get (hv_node, ''Value'’');"); 
snode = 1; | 


hv_ slidtxt = uicontrol (hfv, 'Style','text', 'Position', [20 ZOO 120: 20 }p-se-< 
'String','Time Slider Value: '); 


hv_slidtxt2 = uicontrol (hfv, 'Style', 'text', 'Position',... 
[160 260 80 20], 'String','0O'); 


hv_slidtxt3 = uicontrol (hfv, 'Style', 'text', 'Position',... 
[240 260 80 20],'String',' Seconds"); 


hv_slid = uicontrol (hfv, 'Style', 'slider', 'Position',... 
[20 220 (wd - 40) 20], 'Min',1, 'Max',nsteps, 'Value',1,... 
"Callback',... 
['tstp = round(get(hv_slid, ''Value''));',... 
'set (hv_slidtxt2,''String'',num2str(time(tstp)))7']); 
tstp = 1; 


set (hfv, 'Visib',‘'on"'); 


SESSSELEELELESESSESEEEESES End of CONTROL.M SIIFIFLLFEFFETSSSSSSSSSB 
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% displace.m | : 

% Display the element and its displacement at each time step 

% | McDermott 1999 
SESEESESEESESESSEEEESESSETEEETESESEEEEETEESESEEESSEESEFEEESESSLEEEES 


% plot the element (s) 

if flagd 
zmax=max (z) +max (max (zu))+.1 
ymax=max (y) +max (max (yu) )+.1 
xmax=max (x) +max (max (xu) )+.1; 
zmin=min(z)+min(min(zu))-.1; 
ymin=min (y)+min(min(yu))-.1; 
xmin=min (x) +min (min (xu) )-.1; 
AxMax = [xmin, xmax, ymin, ymax, zmin, zmax] ; 


° 
f 


= 


end; 

m=moviein (nsteps) ; 

 % Display scaling factor in title bar cat 
dispstring = ['Displacement: Scale = ' num2str(dscale, '%1.4g')]}; 

set (gcf, 'Name',dispstring) ; 


for t=l:nsteps 
% display global node numbers 


cla,xlabel('x axis'),ylabel('y axis'),zlabel('Z axis'),grid on,... 


set (gca, 'Box', 'off'),hold on; 


for i=l:ne | 
xl=x (nodes (i,1))+xu(nodes(i,1),t); 
x2=x (nodes (i,2))+xu(nodes(i,2),t); 
x3=x (nodes (1,3) )+xu (nodes (i,3),t); 
x4=x (nodes (i, 4))+xu(nodes(i,4),t); 
yl=y (nodes (i,1))+yu(nodes(i,1),t); 
y2=y (nodes (i,2))+yu(nodes(i,2),t); 
y3=y (nodes (i,3))+yu(nodes(i,3),t); 
y4=y (nodes (1,4) )+yu (nodes (i,4),t); 
zl=z (nodes (i,1))+zu(nodes(i,1),t); 
z2=z (nodes (1,2) )+zu(nodes(i,2),t); 
z3=z (nodes (i, 3))+zu(nodes(i,3),t); 
z4=z (nodes (i, 4))+zu(nodes(i,4),t); 


xvec=[x1,x2,x3,x4,X1]; 


yvec=[yl,y2,y3,y4,yl]; 
zvec=[zl1,z2,2z3,24,21]; 


plot3 (xvec, yvec, zvec) , view (AdjView) , axis (AxMax) ; 
end; 
hold off; 
m(:,t)=getframe; 
end; — 
movie (m) ; 
‘SESSSESSESEESSSSESESSSSSEESE End of DISPLACE.M SSSSEESESESSSEESESEEESSSSS 
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% shownode.m 

S Display the element 

% | McDermott 1999 

SESSELLLLL PSII LISLE LLL EES SESEEEESEEEESESESEESSSSSSFSSESSSSTSESS 


% Remove previous contents 
cla; 


% plot the element (s) 
offset = 0.1 * (max(z) + max(y) + max(x))/3; 


zmax=max (z)+offset; 
ymax=max (y)+offset; 
xmax=max (x)+offset; 
zmin=min(z)-offset; 
ymin=min(y)-offset; 
xmin=min (x) -offset; 
AxMax = [xmin, xmax, ymin, ymax, zmin, zmax]; 


xlabel('x axis'),ylabel('y axis'),zlabel('z axis'),... 
axis (AxMax) , view (AdjView),grid,hoid on; 


for i=l:ne 
xi=x (nodes (i,1)); 
x2=x (nodes (i,2)); 
x3=x (nodes (i,3)); 
x4=x (nodes (i,4)); 
yi=y (nodes(i,1)); 
y2=y (nodes (i,2)); 
y3=y (nodes (i,3)); 
y4=y (nodes (1,4)); 
zl=z (nodes (i,1)); 
z2=z (nodes (i,2)); 
z3=z (nodes (i,3)); 
z4=z (nodes (i,4)); 


xvec=[x1,x2,x3,x4,x1l]; 
yvec=[yl, y2,y3,y4,y1]; 
zvec=[21,z22,23,24,2Z1]; 


plot3 (xvec, yvec, zvec) ; 
end; 


d=offset * 0.25; 
for i=l:nnodes 

text (x(i)td, y(i)+(2*d),z(i)+d,num2str(i)); 
end 


SISSSTESLEETESESELESSESESSEEESESE End Of SHOWNODE.M SESLBSSBBSSSSSSSSSSEB 
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% drawstress. m - plots Von-Mises stress and Yield stress 
% McDermott 1999 
SUGADEGARDEGEATESSERSESRRTTESERESESEREREESERES EER ERESEEEEESEEEEE 


figure (hfa); 
set (hfa,'Visib','on'); 
flaga = 1; 


% tstp = time vector index (not all plots need) 
selected integration eer 

selected element 

% type = Graph Type selected: 


oe 
OS 
O 

b- 

0 

ct 
hot 


% 1 = Stress vs Strain 
% 2 = Stress vs Time 
% 3 = Strain vs Time 
% 4 = Porosity vs Time : 
% 5 = Porosity vs Plastic Strain 
s 6 = Z Stress vs Thickness (at timeé(tstp) ) 
% 7 = X Stress vs Time 
S 8 = Y Stress vs Time 
% 9 = Z Stress vs Time 
S 10 = Stress Multiplot (all 4 stress Siete vs Time) 
% 11 = X Displacement vs Time 
% 12 = Y Displacement vs Time 
% 13 = Z Displacement vs Time 
% 14 = Theta X vs Time 
% 15 = Theta Y vs Time 
S 16 = Theta Z vs Time 
% 17 = Mean Displacement vs Time 
% 18 = Plastic Strain vs Time 
% 19 = Elastic Strain vs Time 
index = (selem-1)*6*ipoint + (point-1)*6; 
indexl = (selem-1)*ipoint + point; 
Switch type 

case l, 

titlestr = ['Element #' num2str(selem) .... 


" Integration Point #' num2str(point)]; 


plot (effstrain(index1,:),vmstress (indexl,:),effstrain(indexl,:),... 


yieldstrs (indexl,:),'--'),grid,xlabel('Effective Strain'),... 


ylabel('Effective Stress'),title(titlestr) ,... 
legend('VME Stress', 'Yield Stress') 

maxvm = max(vmstress (indexl1,:)) 

maxstrn = max(effstrain(indexl,:) ) 


case 2, 


titlestr = ['Element #' num2str(selem) ' Integration Point #', 


num2stxr (point) ]; 

plot (time, vmstress (indexl1,:),time, yieldstrs (indexl]1, :),... 
'-~-'),grid,xlabel('Time (sec) '),... . 
ylabel('Effective Stress'),title(titlestr),... 
legend('VME Stress’, 'Yield Stress"); 

maxvm = max(vmstress (indexl1,:)) 

maxys = max(yieldstrs (indexl1,:)) 
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case 3, 
titlestr = ['Element #' num2str(selem) ' Integration Point #’,.. 


num2str (point) J); 
plot (time, effstrain(index1,:)),grid,xlabel('Time (sec) ") yp... 
ylabel ('Effective Strain'),title(titlestr); 
maxstrn = max(effstrain(indexl, :)) 


case 4, | 
titlestr = ['Element #' num2str(selem) ' Integration Point Ho pak 
num2str (point) J; 
plot (time, void(indexl,:)),grid,xlabel ('Time (sec) '),... 
title(titlestr), ylabel('Porosity'); 
maxpor = max(void(indexl, :)) | ae 


case 5, 
titlestr = ['Element #' num2str(selem) ' Integration Point #',... 
. | num2str (point) ]; 
plot (effplstrn(indexl1,:),void(index1,:)),grid,.. 
xlabel('Effective Plastic Strain'),... 
ylabel('Porosity'),title(titlestr) ; 
maxpor = max(void(indexl]1, :)) 


case 6, 

titlestr = ['Element #' num2str(selem) ' Time: ',.. 
num2str(time(tstp)) " sec']; 

strs = zeros(ipoint,1); 

pnts = l:ipoint; 

for p=l:ipoint 
index3 = (selem-1)*6*ipoint + (p-1)*6 + 3; 
strs(p) = stress(index3,tstp) ; 3 

end 

maxstrs = max(strs(p) ) 

minstrs = min(strs(p)) , 

plot (strs,pnts),grid,axis ij,ylabel('Integration Point’)... 
xlabel('Z Stress'), title(titlestr); 


case 7, 
titlestr = ['Element #' num2str(selem) ' Integration Point #',... 
num2str (point) ]; | 
plot (time, stress (index+1,:)),grid,xlabel('Time (sec) "),.. 
ylabel('X Stress'),title(titlestr); — 
maxstrs = max(stress(indextl,:)) 
minstrs = min(stress (indextl1,:)) 


case 8, | 
titlestr = ['Element #' num2str(selem) ' Integration Point #',... 


num2str (point) ]; 

plot (time, stress (index+2,:)),grid,xlabel ("Time (sec) '),.. 
ylabel('Y Stress'),title(titlestr) ; 

maxstrs = max(stress (indext2,:)) 


minstrs min (stress (indext+2,:)) 
case 9, 
titlestr = ['Element #' num2str(selem) ' Integration Point #',.. 


num2str (point) ]; 
plot (time, stress (index+3,:)),grid,xlabel('Time (sec) ’’ ) yd. 
- ylabel('Z Stress'),title(titlestr); 
maxstrs = max(stress (indext3,:)) 
minstrs = min(stress (indext3,:)) 
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case 10, 


titlestr = ['Element #' num2str(selem) ' Integration Point #',.. 


num2str (point) J; 

subplot (2,2,1) 

plot (time, stress (index+1l,:)),grid,xlabel('Time (sec)'),... 
ylabel('X Stress'),title(titlestr); 

subplot (2,2, 2) 

plot (time, stress (indext+2,:)),grid,xlabel('Time (sec)'),... 
ylabel('Y Stress'), title (titlestr) ; 

subplot (2,2,3) 

plot (time, stress (indext3,:)),grid,xlabel('Time (sec)"),... 
ylabel('Z Stress'),title(titlestr) ; 

subplot (2,2, 4) 

plot (time, vmstress (indexl,:),time, yieldstrs (index1,:),... 
'j--'),grid,xlabel('Time (sec) '),... 
ylabel('Effective Stress'),... 
legend('VME Stress','Yield Stress'),title(titlestr) ; 

maxstrs x = max(stress(index+l, :) ) 

minstrs x = min(stress(index+tl, :)) 

maxstrs y = max(stress(index+2,:)) 

minstrs_y = min(stress (index+2, :)) 

maxstrs 2 = max(stress (indext3, :)) 

minstrs Zz min(stress (index+3,:)) 


case ll, 
titlestr = ['Node #',num2str(snode) ]; 
plot (time, xu(snode,:)/dscale),grid,xlabel('Time (sec) '),... 
ylabel('X Node Displacement'),title(titlestr) ; 
maxdx = max(xu(snode,:))/dscale 
mindx = min(xu(snode,:))/dscale 


case 12, 
titlestr = ['Node #',num2str(snode) ]; 
plot (time, yu(snode,:)/dscale), Sis eet Time (sec)'),... 
ylabel('Y Node Displacement'),title(titlestr); 
maxdy = max(yu(snode,:))/dscale 
mindy = min(yu(snode,:))/dscale 


case 13, 
titlestr = ['Node #',num2str(snode) ]; 
plot (time, zu(snode, :)/dscale),grid,xlabel('Time (sec) '),... 
 ylabel('Z Node Displacement'),title(titlestr) ; 
maxdz = max(zu(snode,:))/dscale 
mindz = min(zu(snode, :))/dscale 


case 14, 
titlestr = ['Node #',num2str(snode) ]; 
plot (time, thxu(snode,:)),grid,xlabel('Time (sec)'), 
ylabel('Theta X'),title(titlestr) ; 
maxdz = max(thxu(snode, :) ) | 
mindz min (thxu (snode, :) ) 


case 15, 
titlestr = ['Node #',num2str(snode) ]; 
plot (time, thyu(snode,:)),grid,xlabel('Time (sec)'),... 
 ylabel('Theta Y'),title(titlestr) ; 
maxdz = max(thyu(snode, :) ) 
mindz = min(thyu(snode, :)) 
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case 16, 
titlestr = ['Node #',num2str(snode) ]; 


plot (time, thzu(snode,:)),grid,xlabel('Time (sec)"), 
ylabel('Theta Z2'),title(titlestr) ; 


maxdz = max(thzu(snode,:) ) 
‘mindz = min(thzu(snode, :) ) 
case 17, 


titlestr = ['Node #',num2str(snode) ]; 

du = sqrt (xu(snode,:).*2 + yu(snode,:).%2 + .. 
zu(snode,:).*2)/dscale; 

maxdu = max (du) 

mindu = min(du) 

plot (time, du) , grid, xlabel ('Time (sec)'),.. 
ylabel('Mean Node Displacement"), tile (eibecens: 


case 18, | 
titlestr = ['Element #' num2str(selem) ' Integration Point #',.. 
num2str (point) ]; 
plot (time, effplstrn(indexl,:)),grid,xlabel('Time (sec)'), 


ylabel ('Effective Plastic Strain'), title (titlestr) ; 
maxplaststrn = max(effplstrn(indexl, :)) 


case 19, 
titlestr = ['Element #', num2str(selem) ' Integration Point #', 
num2str (point) ]; | 
plot (time,estrain(indexl,:)), grid, xlabel('Time (sec)'),... 


ylabel('Effective Elastic Strain' ), title (titlestr); 
maxelaststrn = max(estrain(indexi1,:)) 


otherwise, 
disp('Invalid Graph Type"); 
end 


SESSSSEESEEEETESESESSEESS End of DRAWSTRESS..M &333SSS%SSSSSBSSFBSFBB 
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APPENDIX C. FILE LISTING FOR FEA SHELL ELEMENT SUBROUTINE. 


CR KK RR KK RK RK KK KK KKK KKK KK KERR KER KEKE EKEKEREKKRKEKEEKRKKRKEKKEKEKRKKRKKKRK KKK KK KE 


c* elsh4l1.f 
c* Contains the following subroutines: 
ot elsh4l1 - element matrices for 3-D degenerated 4-node shell 
CR RR RR RK kk ke RR RR KR KR RK RK KR KER KEK IKK RK KR RK KKK RIK KR KK RR 

subroutine elsh4l (nnds,xcoord, ycoord, zcoord,nesh4,ndsh4,matsh4, 

1 nmat,matno, eprop, ndsdof, sdisp, svel, sforce, 

2 smass,nsdof,ishel4,stssh4,mpoint, ipoint,itime) 
CER RKKK KKK KKK KKK KK KKK KKK KKK KKKEKKKKKKKKRKKKKKKKKRAEKKKKKKKKREK KK KKK KKK KKK KE 
compute and assemble element matrices for shell elements with 
4 nodes (linear analysis) 
3-D degenerated formulation including transverse shear and. 
transverse normal deformation 
pressure formulation 


qaqaqgnana 


McDermott 1999 
CEKKE KKK KER KKK KKK KKK KKK KKK KKK KEKE KKK KKK KKK KKK KEKE KEKE KKK RE RK KKK KR KK 
common /fblk/ fail,arrayl1 (125) 
dimension xcoord(nnds), ycoord(nnds),zcoord(nnds), 
ndsh4 (nesh4+1,4),matsh4 (nesh4+1),ishel4 (nesh4+1), 
matno (nmat),eprop (nmat, 40),ndsdof(nnds), 
sdisp (nsdof),svel (nsdof),sforce(nsdof),smass(nsdof), 
stssh4 (mpoint,nesh4*8) 
dimension eforce(24),lnode(4),evel(24),slope(10),dstrnplist (6), 


& WD 


1 emass (24),gausx(4),lindex(24),strain(10),corr(2), 
2 array2 (15),estrain(6),estress(6),s(6),A(2,2),Ainv(2,2), 
3 gausy (4) ,weighx (4) ,weighy(4),shapel(3),derivl1(3), 
4 shapef (4) ,derivl(2,4),derivg(3,4),aj(3,3),ajinv(3,3), 
5 bmtx (6,24) ,tmp(8),xc(4),yc(4),zc(4),estressg(6), 
6 derilg(3,4),d1(3),edisp(24),edispg(24),rot6t(6,6), 
7 d2(3),v1(3),v2(3),v3(3),gausz(10),weighz(10), 
8 bmtxt (24,6),delt(6),t1(3,3),tlinv(3,3),B(2), 
9 rot6 (6,6) ,rot3(3,3),estrainp(6),eforceg(24), 
1 weighxh (4) ,weighyh (4) ,weighzh(4),eforceh(24), 
2 gausxh (4),gausyh(4),gauszh(4),iplane(5) 
Cc : 
data ngausx,ngausy,ndof,n /1, 1, 6, 4/ 
data ngausxh,ngausyh,ngauszh,hourl /2, 2, 1, 0.05/ 
data ref,eps /0.0,1.0e-4 / 
data delt /1.0, 1.0, 1.0, 0.0, 0.0, 0.0/ 
data iotene Wl, 2; Ag Sy. G/ 
neldof=ndof*n 
ngausz=ipoint 
ngauss=ngausx*ngausy*ngausz 
ngaussh=ngausxh*ngausyh*ngauszh 
S 


c *** extract gauss points and weights for numerical integration 


call gauss (ngausx, gausx, weighx) 
call gauss (ngausy, gausy, weighy) 
call gauss (ngausz, gausz,weighz) 
call gauss (ngausxh, gausxh, weighxh) 
call gauss (ngausyh, gausyh, weighyh) 
call gauss (ngauszh, gauszh,weighzh) 
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‘i | 
c *** Joop for number of elements 
Cc 
do 370 ie=1,nesh4 
CS . 
igausp=0 
eS 
c *** initialize element internal force vector 
c 
do 10 i=1,neldof 
eforceh(i) = 0.0 
eforce(i) = 0.0 
10 continue 
c 
c *** extract material property index and material properties 
és | 
imt=matsh4 (ie) 
mattyp=matno (imt) 
densty=eprop(imt,1) . 
elast = eprop(imt,2) 
poiss = eprop(imt, 3) 
thick = eprop(imt, 4) 
Syp = eprop(imt, 5) 
ql = eprop(imt, 6) 
q2 = eprop(imt,7) 
q3 = eprop(imt, 8) 
fn = eprop(imt, 9) 
en = eprop(imt,10) 
sn = eprop(imt,11) 
PhiO = eprop(imt,12) 
iseg = eprop(imt,13) + 0.1 
ystrain = Syp / elast 
strain(1) = ystrain 
if (iseg.gt.0) then 
do 20 i=1,iseg 
slope(i) = eprop(imt,12 + i * 2) 
strain(itl) = eprop(imt,13 + i * 2) 
20 continue 
endif 
bulk = elast / (3.0 * (1.0 - 2.0 * poiss)) 
shear = elast / (2.0 * (1.0 + poiss)) 
C 
c *** Default Values 
c 
if (ql .ne. 0.0) then 
if (q2 .eq. 0.0) q2 = 1.0 
if (q3 .eq. 0.0) g3 = qli**2 
if (fn .ne. 0.0) then 
if (sn .eg. 0.0) sn = 0.1 
if (en .eq. 0.0) en = 0.3 
endif : 
endif 
c 
c *** find dof index of each element 
G 


do 30 i=1,n 
| lnode (i) =ndsh4 (ie, i) 
do 30 j=1,ndof 
lindex ((i-1) *ndof+)j) 


30 continue 
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=ndsdof (lnode(i))+}-1 
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load displacement for element 


do 40 i=1,neldof 

edispg(i) = sdisp(lindex(i)) 
evel(i) = svel (lindex({i)) 
continue 


extract element nodal coordinate values 


do 50 i=l1,n 

xc (1) =xcoord (lnode (i) ) 
yc (1) =ycoord (lnode (1) ) 
zc (i)=zcoord (lnode (i) ) 
continue 


compute direction vectors 
(vl & v2: tangent, v3: normal to surface) 


ai (1) =xe(2)-xce (1) 
dl (2) =ye (2) -ye(1) 
dl (3) =zc(2)-ze (1) 
d2 (1) =xe (4) -xc (1) 
a2 (2) =ye (4) -yc(1) 
d2 (3) =ze(4)-zc(1) 


Calculation of Hourglass Control Multiplier 


Need Surface Area of Element 
all=sqrt (d1 (1) **2+d1 (2) **2+di (3) **2) 


al4 = sqrt (d2(1)**2+d2 (2) **2+d2 (3) **2) 
al2 = sqrt ((xce(3)-xe(2) ) **2+(ye(3)-ye (2) ) **2+(zc(3)-zce(2) ) **2) 
al3 = sqrt((xe(3)-xe (4) ) **2+ (yc(3) -~yo (4) ) **2+(ze(3)-ze(4))**2) 


adotl = (xc(1)-xc(2))*(xe(3)-xe(2))+(yce(1)-ye(2) )* (ye(3)-ye (2) ) 
+ (zc0(1)-zc(2))*(ze(3)-zce(2)) 

adot2 = (xc(1)-xec(4))*(xc(3)-xe(4))+(ye(1) ~ye(4))* (ye (3) ~yeo(4)) 
+ (ze(4) -z0c(4))*(zce(3)-zc(4) ) 

area = 0.5*(sgrt(all**2*al2**2-adot1**2) + 

sart (al3**2*al4**2-adot2**2) ) 

hour = hourl * thick**2 / area 


v3 is normal to plane 
v3 (1) =d1 (2) *d2 (3) -d1 (3) *d2 (2) 
v3 (2) =d1 (3) *d2 (1) -d1 (1) *d2 (3) 
v3 (3) =dl1 (1) *d2 (2) -d1 (2) *d2 (1) 
sum=sqrt (v3 (1) **2+v3 (2) **2+v3 (3) **2) 
v3(1)=v3(1)/sum 
v3 (2)=v3 (2) /sum 
v3 (3)=v3 (3) /sum 


v2 is unit vector in direction of node 1 to node 4 
~w2(1) = d2(1)/al4 
v2(2) = d2(2)/al4 
v2(3) = d2(3)/al4 


vl is orthogonal to v2, in plane normal to v3 
vl (1) =v3 (2) *v2 (3) -v3 (3) *v2 (2) 
vl (2) =v3 (3) *v2 (1) -v3 (1) *v2 (3) 
vl (3) =v3 (1) *v2 (2) -v3 (2) *v2 (1) 
sum = sqrt (vl(1)**2+vl (2) **2+v1 (3) **2) 
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Cc 


vl (1)=-v1(1)/sum 
vl (2)=-v1(2)/sum 
v1(3)=-v1 (3) /sum 


c *** Strain Transformation Matrix 


Q 


qaa 0 


Q 


kkk 


60 


kkk 


70 


wk eK 


Rotation tranformation matrix, 


rot6é(1,1)=v1(1)**2 

roct6é (1,2) =v1 (2) **2 
rot6(1,3)=vl1(3)**2 
rot6é(1,4)=vl1 (1) *vl (2) 
rot6(1,5)=v1 (2) *v1 (3) 

roté6é (1, 6)=v1 (1) *vl1 (3) 
rot6(2,1)=v2 (1) **2 

rot6 (2,2) =v2 (2) **2 

rot6 (2,3) =v2 (3) **2. 
rot6(2,4)=v2 (1) *v2 (2) 
rot6(2,5)=v2 (2) *v2 (3) 

rot6 (2,6) =v2 (1) *v2 (3) 

roté6 (3,1)=v3(1)**2 

rot6é (3,2)=v3 (2) **2 

roté6 (3,3) =v3 (3) **2 

rot6 (3, 4)=v3 (1) *v3 (2) 
rot6(3,5)=v3 (2) *v3 (3) 

rot6 (3, 6)=v3 (1) *v3 (3) 
rot6(4,1)=2.0*vi (1) *v2(1) 
rot6(4,2)=2.0*vl1 (2) *v2 (2) 

rot6 (4,3)=2.0*vil (3) *v2 (3) 
rot6(4,4)=v1 (1) *v2(2)+v2(1) *v1(2) 
roté6é (4,5) =vl1 (2) *v2 (3) +v2 (2) *v1 (3) 
rot6 (4,6) =v1 (3) *v2 (1) +v2 (3) *v1 (1) 
rot6(5,1)=2.0*v2 (1) *v3(1) 
rot6(5,2)=2.0*v2 (2) *v3 (2) 

roté (5,3)=2.0*v2 (3) *v3 (3) 


—rot6é (5,4) =v2 (1) *v3 (2) +v3 (1) *v2 (2) 


rot6(5,5)=v2 (2) *v3 (3) +v3 (2) *v2 (3) 
roté (5, 6) =v2 (3) *v3 (1) +v3 (3) *v2 (1) 


 roté6(6,1)=2.0*v3 (1) *vi (1) 


rot6 (6,2)=2.0*v3 (2) *v1(2) 

rot6 (6,3)=2.0*v3 (3) *vi (3) | 
rot6 (6, 4)=v3 (1) *v1(2)4+vl (1) *v3 (2) 
rot6 (6,5)=v3 (2) *v1 (3)+vl1 (2) *v3 (3) 
roté (6, 6) =v3 (3) *v1 (1) +v1 (3) *v3 (1) 


Transpose of Transformation Matrix 


do 60 i=1,6 
do 60 j=l, 6 

rotét (i,j) =rot6 (j,i) 
continue 


do 70 i=1,3 
ti(i,1) = vl1(i) 
t1(i,2) = v2(i) 
t1(i,3) = v3(i) 
continue 


call inv3x3(t1,tlinv,det) 


Transform Rotation displacement to local coordinates 
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and its inverse 





do 80 i=1,4 
do 80 j=1,3 
edisp ( (i-1) *6+3j)=edispg ( (i-1) *6+3) 
edisp ( (i-1) *6+3+j) =edispg ( (i-1) *6+4) *tlinv(j,1) 
1 +edispg ( (i-1) *6+5) *tlinv(4j,2) 
2 _  +edispg( (i-1) *6+6) *tlinv(j,3) 
80 _ continue 





c 
c *** Begining of Hourglass control 
fe 
if (hour.gt.0.0) then 
do 140 ix=1,ngausxh 
re=gausxh (ix) 
wx=weighxh (ix) 
do 140 iy=1,ngausyh 
sc=gausyh (iy) 
wy=weighyh (iy) 

do 140 iz=1,ngauszh 

tc=gauszh (iz) 

wz=weighzh (iz) 
c 
c *** compute shape functions and derivatives 
! 

| call shp2d4 (rc,sc,shapef) 

call der2d4 (rc,sc,derivl) 

call shpld2 (tc, shapel) 

call derld2 (tc,deriv1l1) 
C = 
*** compute jacobian matrix and its inverse 


Q 


hz=thick*0.5* (shapel (2) * (1.0-ref)-shapel (1)*(1.0+ref) ) 
hzdt=thick*0.5 


aj (1,1)=derivl (1,1) *xc(1)+derivl (1,2) *xc(2)+derivl (1,3) *xc(3)+ 
Al derivl (1,4) *xc(4)+derivl (1,1) *hz*v3(1)+ 

2 derivl (1,2) *hz*v3(1)+derivl (1,3) *hz*v3(1)+ 

3 derivl (1,4) *hz*v3(1) 


aj (2,1)=derivl (2,1) *xc(1)+derivl (2,2) *xc(2)+derivl (2,3) *xc(3)+ 
Z | derivl (2,4) *xc(4)+derivl (2,1) *hz*v3(1)+ 
2 derivl (2,2) *hz*v3(1)+derivl (2,3) *hz*v3(1)+ 
3 derivl (2,4) *hz*v3(1) 
aj (3,1)=hzdt*v3 (1) 


aj (1,2)=derivl (1,1) *yc(1)+derivl (1,2) *yc(2)+derivl (1,3) *yc(3)+ 
1 derivl (1,4) *yc(4)+derivl (1,1) *hz*v3(2)+ 

2 derivl (1,2) *hz*v3(2)+derivl (1,3) *hz*v3(2)+ 

3 derivl (1,4) *hz*v3(2) 


aj (2,2) =derivl (2,1) *yc(1)+derivl (2,2) *yc(2) +derivl (2,3) *yc(3)+ 
1 derivl (2,4) *yc(4)+derivl (2,1) *hz*v3(2)+ 
2 ) derivl (2,2) *hz*v3 (2)+derivl (2,3) *hz*v3 (2) + 
S derivl (2,4) *hz*v3 (2) 
aj (3,2) =hzdt*v3 (2) 


aj (1,3) =derivl (1,1) *ze(1)+derivl (1,2) *zc(2) +derivl (1,3) *zc(3)+ 
1 derivl(1,4)*zc(4)+derivl (1,1) *hz*v3(3)+ 

2 derivl (1,2) *hz*v3 (3) +derivl (1,3) *hz*v3(3)+ 

3 derivl (1, 4) *hz*v3 (3) 
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Ee eee 


Cc 


Cc 
Cc 
Cc 


I a ae el a Ache a cra re scar a 
1 derivl (2,4)*zc(4)+derivl (2,1) *hz*v3(3)+ 
2 derivl (2,2) *hz*v3 (3)+derivl (2, eae Nake 
3 derivl (2,4) *hz*v3 (3) 
aj (3,3) =hzdt*v3 (3) 


call inv3x3 (aj,ajinv,det) 
*** compute global derivatives and strain-nodal displacement matrix 
do 90 i=1,n | 


derivg(1,i)=ajinv(1, 1)*derivl(1,i)+ajinv(1,2)*derivl (2,1) 
derivg(2,i)=ajinv(2,1)*derivl(1,i)+ajinv(2,2)*derivl (2,i) 
derivg(3, i) =ajinv (3, 1) *derivl(1,i)+ajinv (3,2) *derivl (2,1) 
derilg(1,i)=ajinv(1,3) *hzdt 
derilg(2,1)=ajinv (2,3) *hzdt 
derilg(3, aa 3) *hzdt 
90 continue 


do 100 i=1,6 
do 100 j=1,24 
bmtx (i,3)=0.0 
100 continue 


do 110 i=i,n 
il=(i-1) *6+1 
12=il1tl 
13=12+1 
i4=13+1 
15=i44+1 
i16=i5+1 
gkl=derivg(1,i) *hz+shapef (i) *derilg(1,i) 
gk2=derivg(2,i)*hz+shapef (i) *derilg(2,i) 
gk3=derivg(3,i) *hz+shapef (i) *derilg(3,i) 


bmtx(1,i1)=derivg(1i,i) 
bmtx (1,14)=gk1* (-v2(1)) 
bmtx (1,15)=gki*vi(1) 
bmtx (1,16) =gki*v3 (1) 
bmtx(2,i2)=derivg(2,i) 
bmtx (2,14) =gk2* (-v2 (2) ) 
bmtx (2,15) =gk2*vi1 (2) 
bmtx (2,16) =gk2*v3 (2) 
bmtx (3,13) =derivg (3,1) 
bmtx (3,14)=gk3* (-v2 (3) ) 
bmtx (3,15) =gk3*vi1 (3) 
bmtx (3,16) =gk3*v3 (3) 
bmtx (4,11)=derivg(2,i) 
bmtx (4,1i2)=derivg(i,i) 
bmtx (4, i4)=gk2* (~v2 (1) ) +gk1* (-v2(2)) 
bmtx (4,15) =gk2*v1 (1) +gki*vi (2) 
bmtx (4,16) =gk2*v3 (1) +gk1*v3 (2) 
bmtx (5,1i2)=derivg(3,i) 
bmtx (5,13) =derivg(2,i) 
bmtx (5,14) =gk3* (-v2 (2) ) +gk2* (-v2 (3) ) 
bmtx (5,15) =gk3*vl1 (2) +gk2*vi1 (3) 
bmtx (5,16) =gk3*v3 (2) +gk2*v3 (3) 
bmtx (6,i1)=derivg(3,i) 
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bmtx (6,13)=derivg(1,i) 

bmtx (6,14) =gk3* (~v2 (1) )+gk1* (-v2 (3) ) 

bmtx (6,15) =gk3*vil (1) +gk1*vi1 (3) 

bmtx (6,16) =gk3*v3 (1) +gk1*v3 (3) 
continue 


do 120 i=1,6 
do 120 j=1,neldof 
bmt xt (j,1i)=bmtx (i,j) 


weight = wx * wy * wz 
detwt = det * weight 


calculate strain 

call mmult(6,24,1,bmtx, edisp, estrainp,0,1.0) 
Transform strain to local coordinates System - 

call MAOTOUE NG A, Rote xsteatnn, 6sevain, 05.0) 
Subtract any prior plastic strain 


do 130 i=1,6 
estrain(i) = estrain(i) - stssh4(12+1i,issts) 


Calculate stress using plane-strain formulas 


elast/(1.0-poiss**2)*(estrain(1) + poiss*estrain(2) ) 
elast/(1.0-poiss**2)*(poiss*estrain(1) + estrain(2) ) 
estress (3) elast * estrain (3) 
estress (4) shear * estrain(4) 
- estress (5) 5.0 /6.0 * shear * estrain(5) 
estress (6) 5.0 /6.0 * shear * estrain(6) 


estress (1) 
estress (2) 


Transform stress back to global coordinates 
call mmult (6, 6,1,roté6t, estress, estressg,0,1.0) 

find eforceh for this gauss point and sum to element force 
call mmult (24,6,1,bmtxt, estressg, eforceh, 1,detwt) 


continue 
endif 


end of hourglass loop 


numerical integration loop for membrane 


icount=0 
do 300 ix=1,ngausx 
rc=gausx (ix) 
wx=weighx (ix) 
do 300 iy=1,ngausy 
sc=gausy (iy) 
wy=weighy (iy) 
do 300 iz=1,ngausz 
tc=gausz (iz) 
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wz=weighz (iz) 


increase counter 


igausp=igausp+1 
issts=(ie-1)*8+igausp 


Initialize stssh4 if first time 


if (itime . 


eq. 1) then 


do 150 i=13,21 


stssh4(19,issts) 
stssh4(20,issts) 


endif 


stssh4(i,issts) = 0.0 


SYP 
Phid 


Hh fl 


compute shape functions and derivatives 


call 
call 
call 
call 


compute jacobian 


hz=thick*0.5* (shapel (2) * (1.0-ref)-shapel (1) *(1.O0+ref) ) 


shp2d4 (rc,sc,shapef) 
der2d4 (rc,sc,derivl) 
shpld2 (tc,shapel) 
Gerld2 (tc,derivl) 


matrix and its inverse 


hzdt=thick*0.5 


aj (3, 


aj (3, 


aj (1,1)=derivl (1,1) *xe(1)+derivl (1,2) *xc(2)+derivl (1,3) *xe(3)+ 


derivl (1,4) *xc(4)+derivl (1,1) *hz*v3(1)+ 
derivl (1,2) *hz*v3 (1)+derivl (1,3) *hz*v3(1)+ 
derivl (1,4) *hz*v3(1)_ 


aj (2,1)=derivl (2,1) *xc(1)+derivl (2,2) *xc(2)+derivl (2,3) *xc(3)+ 


derivl (2,4) *xec(4)+derivl (2,1) *hz*v3(1)+ 
derivl (2,2) *hz*v3(1)+derivl (2,3) *hz*v3(1)+ 
derivl (2,4) *hz*v3 (1) 

1)=hzdt*v3 (1) 


aj (1,2)=derivl(1,1)*yc(1)+derivl (1,2) *yc (2) tderivl (1,3) *yc(3)+ 


derivl (1,4) *yc(4)+derivl (1,1) *hz*v3(2)+ 
derivi (1,2) *hz*v3 (2)+derivl (1,3) *hz*v3(2)+ 
derivl (1, 4) *nhz*v3 (2) : 


aj (2,2)=derivl (2,1) *yc(1)+derivl (2,2) *yc(2)+derivl (2,3) *yc(3)+ 


derivl (2,4) *yc (4) +derivl (2,1) *hz*v3(2)+ 
derivl (2,2) *hz*v3(2)+derivl (2,3) *hz*v3(2)+ 
derivl (2,4) *hz*v3 (2) 
2) =hzdt*v3 (2) 


aj (1,3)=derivl (1,1) *zc(1)+derivl (1,2) *zc(2)+derivl (1,3) *zce(3)+ 


derivl (1, 4)*zc(4)+derivl (1,1) *hz*v3(3)+ 
derivl (1,2) *hz*v3(3)+derivl (1,3) *hz*v3(3)+ 
derivi (1,4) *hz*v3 (3) 


aj (2,3) =derivl (2,1) *zc(1)+derivl (2,2)*ze(2)+derivl (2,3) *zce(3) +. 


derivl (2,4)*zc(4)+derivl (2,1) *hz*v3 (3) + 
derivil (2,2) *hz*v3 (3) +derivl (2,3) *hz*v3(3)+ 
derivl (2,4) *hz*v3 (3) 
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aj (3,3) =hzdt*v3 (3) 
call inv3x3 (aj,ajinv, det) 

*** compute global derivatives and strain-nodal displacement matrix 
do 160 i=i,n 


derivg(1,i)=ajinv(1,1)*derivl (1,1) +ajinv(1,2) *derivl (2,i) 
derivg(2,i)=ajinv(2,1)*derivl (1,i)+ajinv (2,2) *derivl (2,i) 
derivg (3,i)=ajinv(3,1)*derivl(1,i)tajinv (3,2) *derivl (2,1) 
derilg(1,i)=ajinv (1,3) *hzdt 
derilg(2,i)=ajinv(2,3)*hzdt 
derilg(3,1i)=ajinv (3,3) *hzdt 
160 continue 


do 170 i=1,6 
do 170 j=1,24 
bmtx (i, 3) =0.0 
170 continue 


do 180 i=1,n 

i11=(i-1) *6+1 

i2=il1+1 

13=12+1 

14=13+1 

i5=14+1 

16=15+1 
gki=derivg(1,i)*hz+shapef(i)*derilg(1,i) 
gk2=derivg(2,i) *hz+shapef (i) *derilg(2,i) 
gk3=derivg(3,i) *hz+shapef (i) *derilg(3;,21) 


bmtx(1,i1)=derivg(1,i) 
bmtx (1,14) =gk1* (-v2 (1) ) 
bmtx (1,i15)=gk1*vi (1) 
bmtx (1,16) =gk1*v3 (1) 
bmtx (2,i2)=derivg(2,i) 
bmtx (2,14) =gk2* (-v2 (2) ) 
bmtx (2,15) =gk2*v1 (2) 
bmtx (2,16) =gk2*v3 (2) 
bmtx (3,13) =derivg (3,1) 
bmtx (3,14) =gk3* (-v2 (3) ) 
bmtx (3,15)=gk3*vi (3) 
bmtx (3,16) =gk3*v3 (3) 
bmtx (4,i1)=derivg(2,i) 
bmtx (4,i2)=derivg(1,i) . 
bmtx (4,14) =gk2* (-v2 (1) )+qk1* (-v2 (2) ) 
bmtx (4,15) =gk2*vi1 (1) +gk1*vl (2) 
bmtx (4,16) =gk2*v3 (1) +gk1*v3 (2) 
bmtx (5,12)=derivg (3,1) 
bmtx (5,13)=derivg(2,i) 
bmtx (5,14) =gk3* (-v2 (2) )+gk2* (-v2 (3) ) 
bmtx (5,15) =gk3*vl1 (2) +gk2*v1 (3) 
bmtx (5,16) =gk3*v3 (2) +gk2*v3 (3) 
bmtx (6,i1)=derivg (3,1) 
bmtx (6,13) =derivg(1,i) 
bmtx (6,14) =gk3* (~v2 (1) ) +gk1* (-v2.(3) ) 
bmtx (6,15)=gk3*v1(1)+gk1*vl1 (3) 
bmtx (6,16) =gk3*v3 (1) +gk1*v3 (3) 
180 continue 
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*** calculate transpose of bmtx 


do 190 i=1,6 
do 190 j=1,neldof 
190 ' bmtxt (j,i)=bmtx (i,j) 


weight = wx * wy * wz 
detwt = det * weight 
icount=icountt+1l 
tmp (icount) =detwt 


*** calculate strain 

call mmult (6,24,1,bmtx, edisp, estrainp, 0,1.0) 
*** Transform strain to local coordinates system 

call cnaetGne aneeeb eeeea un estenin.0/4:.0) 
ae Subtract any prior plastic strain 


do 200 i=1,6 
200 estrain(i) = estrain(i) - stssh4(12+i,issts) 


stressO = stssh4(19,issts) 
Phi = stssh4(20,issts) 
plastrn = stssh4(21,issts) 


*** Calculate stress using plane-strain formulas 


= elast/(1.0-poiss**2)*(estrain(1) + poiss*estrain (2) ) 
= elast/(1.0-poiss**2)*(poiss*estrain(1) + estrain(2)) 
estress (3) elast * estrain(3) 

estress (4) shear * estrain(4) 

estress (5) 5.0 /6.0 * shear * estrain(5) 
estress (6) 5.0 /6.0 * shear * estrain(6) 


estress(1) 
estress (2) 


how tt tl 


*** Elastic-Plastic Calculations 


p = -(estress(1) + estress(2) + estress(3)) / 3.0 
do 210 i=1,6 | 
210 s{i) = estress(i) + p * delt(i) 
q = sqrt(1.5 *(s(1)**2 + s(2)**2 + s(3)**2 + 2.0 * s(4)**2 + 
i 2.0 * s(5)**2 + 2.0 * s(6)**2)) 
F = (q / stress0) **2+2.0*q1l*Phi*cosh (-1.5*q2*p/stress0) - 
1 (1.0+q3* Phi**2) 
if (F.gt.Q0) then 
iter = 0 
dep = 0.0 
deq = 0.0 
** Begin Iteration Loop ** 
220 cl = 3.0*q2/(2.0*stress0**2) 


Fql = 2.0 * q / (stress0**2) 
Fq2 = 2.0 / (stress0**2) 
100 








230 


Fpl = 2.0 * qi * Phi * -cl * sinh(-cl * p) 
Epa = 250. G1. Pri. * el**2 * cosh (=-cl *p) 
Fsl = -q **2 / (2.0 * stress0**3) + 2.0%*g1*Phi*cl* 


p*sinh(-cl*p) 
Fsp = 2.0*q1l*Phi* (-cl**2 * p*cosh(- ae + 
cl/stressO * sinh(-cl*p)) 
Fsq = -q / stress0**3 


stressl = Syp 
esl = stress0/elast | 
el = -p*dep/((1-Phi) *stress0) 
psl = plastrn + eil + esl 
iflag = 0 
do 230 i=1l,iseg 

if (iflag .eq. 0) then 

if(psi .it. strain(it1l)) then 

stressl = stressi + slope(i)*(psl-strain(i) ) 


= i 
iflag = 1 
stressl = stressl + slope(i)*(strain(itl)-strain(i)) 
endif i 


endif 
continue 
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c 
c *** Numerical Error Trap 
C 
if ((stressl .eq. stress0) .or. (abs(el) .eq. 0.0)) 
then 
dsdep = slope(isegl) 
else 
dsdep = (stressl-stress0) /el 
endif 
if (abs(dsdep) .gt. elast) then 
| dsdep = slope(isegl) 
endif 
Cc 


stressl = Syp 
el = q * deq / ((1-Phi) * stress0) 
psl = plastrn + el 
iflag = 0 
do 240 i=1,iseg 

if (iflag .eq. 0) then 

if(psl .1t. strain(itl)) then 

stressl = stressl + slope(i)*(psi-strain (i) ) 


isegl = i 
iflag = 1 
else 
stressl = stressl + slope(i)*(strain(itl)- 
strain (i) ) “s 
endif 
endif. 
240 continue 
c 
c *** Numerical Error Trap 
c 
if ((stressl .eq. stress0).or.(abs(e1l) .eq. 0.0)) then 
dsdeq = slope(isegl) 
else 
dsdeg = (stressl-stress0)/el 
endif 
if (abs(dsdeq) .gt. elast) then 
dsdeq = slope(isegl) 
endif 
c 
c *** Using only stress0O derivatives 
‘ | | 
A(1,1) = bulk * Fpl1+Fsl*dsdep 
A(1,2) = -3.0*shear*Fql+Fsl*dsdeq 
A(2,1) = Fql+deq* (bulk*Fp2+Fsp*dsdep) +dep*Fsq*dsdep 
A(2,2) = Fpl+dep* (-3.0*shear*Fq2+Fsq*dsdeq) +deq*Fsp*dsdeq 
Cc 
B(1l) = -1.0 * F 
B(2) = -1.0 * dep * Fql - deq * Fpl 
c 
call inv2x2 (A, Ainv,d) 
Cc : 
call mmult(2,2,1,Ainv,B,corr,0,1.0) 
c 
Phi = stssh4(20,issts) 
plastrn = stssh4(21,issts) 
stressO = stssh4(19,issts) 
Cc 


dep = dep + corr(1) 
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Q 





deq = deq + corr(2) 


degrowth = 0.0 

do 250 i=1,5 
ind = iplane(i) 
dstrnplst (ind) = dep * delt(ind) / 3.0 + 
deq *°3.0 * s(ind) / (2.0 * q) 

continue 


dstrnplst (3) = dep / 3.0 + deq * estress(3) / q 


do 260 i=1,3 
260 degrowth = degrowth + dstrnplst (i) 


** Calculate New stress, p, q, and F 


estress(1) = elast/(1. 0-poiss**2) *(estrain(1) =i deeeisieuiii + 
- | poiss* (estrain(2)-dstrnplst (2) )) 
, estress(2) = elast/(1.0-poiss**2) * 
al (poiss* (estrain(1)-dstrnplst (1) )+estrain (2) ~dstrnplst (2) ) 
estress(3) = elast * (estrain(3)-dstrnplst (3) ) 


estress(4) = shear * (estrain(4)-dstrnplst (4) ) 
estress(5) = 5.0 /6.0 * shear * (estrain(5)-dstrnplst (5) ) 
estress(6) = 5.0 /6.0 * shear * (estrain(6)-dstrnplst (6) ) 

p = (estress(1) + estress(2) + estress(3)) / -3.0 


do 270 i=1,6 
270 s(i) = estress(i) + p * delt(i) 


deltep = (q*deq - p * dep) /((1-Phi) *stress0) 


Prevent reduction in plastic strain 
if (deltep .1t. 0.0) then 
deltep = 0.0 | 
endif 


plastrn = plastrn + deltep 
if (fn .ne. 0.0) then 
anucl = fn/(sn*2.50662827463) *exp (-.5*( (plastrn-en) /sn) **2) 
else 
anucl = 0.0 
endif 


Phil = Phi + (1 - Phi) * degrowth + anucl * deltep 


Phi cannot decrease 


af (PRI: «gt... Phi) -Phi Phil 


Turn off Void Effects if ql = 0.0 
if (ql .eq. 0.0) Phi = 0.0 


Calculate Strain Hardening 
elastrn = stress0/elast 
stressO = Syp 
iflag =0 . 
do 280 i=1,iseg 
if (iflag .eq. 0) then 
psl = plastrn + elastrn 
1f(psl .1lt. strain(itl)) then 
stressO = stress0O + slope (i)*(psl- strain (i)) 
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280 


Kk* 


kk * 


Kk* 


k*k* 


290 


kk* 


300 


k*KK 


310 


w*k* 





iflag = 1 
else | 
stressO = stress0O + slope(i)*(strain(it+t1l)-strain(i) ) 
endif 
endif 
continue 


sqrt(1.5 *(s(1)**2 + s(2)**2 + s(3)**2 + 2.0 * s(4)**2 + 
1 | 2.0 * s(5)**2 + 2.0 * s(6)**2)) 


1Q 
ll 


(q / stress0) **2+2.0*q1*Phi*cosh (-1.5*q2*p/stress0) - 
i (1.0+q3* Phi**2) 


ty} 
rl 


iter = iter +l 

if (iter.gt.500) then 
write(*,*) 'Failed to Converge’ 
stop | 


endif 
if(abs(F).gt.eps) goto 220 


endif 
End of Plastic Deformation Section 
Transform stress back to global coordinates 

call mmult (6,6,1,rot6t,estress, estressg, 0,1.0) 
find eforce for this gauss point and sum to element force 
| call mmult (24,6,1,bmtxt, estressg, eforce, 1, detwt) 
save stress and strain 


do 290 i= 1,6 
stssh4(i,issts)=estress (i) 


stssh4 (it+6,issts)=estrain(i) - dstrnplst (i) 
stssh4(i+12,issts)=stssh4(i+12,issts) + dstrnplst (i) 
continue 
stssh4(19,issts) = stress0 
stssh4(20,issts) = Phi 
stssh4(21,issts) = plastrn 


end of integration loop 
continue 
Hourglass Control Correction 


do 310 i= 1,24 
eforce(i) = eforce(i) + hour * (eforceh(i) - eforce(i)) 


Transform Moments back to Global Coordinate system 


eforceg((i-1) *6+j) =eforce ((i-1) *6+j) 
eforceg ( (i-1) *6+3+3) =eforce ((i-1) *6+4) *t1(j,1) 
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Q 


@) 


re) 


kK 


340 


350 


kk* 


360 


kkk 


370 


comput 


ems=s 


assemb 


end of 





+eforce ( (i-1) *6+5) *t1 (4,2) 
+eforce ( (i-1) *6+6) *t1 (4,3) 
continue 


e lumped mass matrix 
sum=0.0 

do 340 i=l, ngauss 
sum=sumt+tmp (i) 


um/n 


do 350 i=l, neldof 
emass (1)=ems*densty 


le mass and force matrices into system matrix 

do 360 i=1, neldof 

sforce(lindex(i)) = sforce(lindex(i)) + eforceg(i) 
smass(lindex(i)) = smass(lindex(i)) + emass(i) 
continue 


element loop 


continue 


return 


end 
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